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TITLE OF THE INVENTION 



SYSTEM AND METHOD OF WAVEFRONT SENSING 



CROSS-REFERENCE TO RELATED APPLICATIONS 



This application claims priority under 35 U.S.C. § 1 19 to co-pending U.S. Provisional 
Patent Application No. 60/394,232, filed on 9 July 2002 in the name of Daniel M. Topa, the 
entirety of which is hereby incorporated by reference herein for all purposes as if fiiUy set forth 
herein. 



BACKGROUND AND SUMMARY OF THE INVENTION 



[0001] Technical Field. 

[0002] This invention pertains to the field of wavefront sensing methods and devices, and more 
particularly, to a method and system for more precisely determiiiing the location of a focal spot in 
a wavefront sensor, for more precisely determining the location of vertices or boundaries of 
lenslets in a lenslet array with respect to pixels of a wavefront sensor, for using this information 
to reconstruct an arbitrary wavefront with the wavefront sensor, and for more efficiently 
performing the calculations to perform the reconstruction. 
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[0003] Description. 

[0004] A light wavefiront may be defined as the virtual surface delimited by all possible rays 
having an equal optical path length from a spatially coherent source. For example, the wavefront 
of light emanating from a point Ught source is a sphere (or a partial sphere where Ught from the 
point light source is limited to emission along a small range of angles). Meanwhile, the 
wavefront created by a coUimating lens mounted at a location one focal length away from a point 
source is a plane. A wavefront may be planar, spherical, or have some arbitrary shape dictated by 
other elements of an optical system through which the light is transmitted or reflected. 
[0005] A nimiber of different wavefront sensors and associated methods are known. Among 
these are the Shack-Hartmann wavefront sensor, which will now be described in greater detail. 
[0006] FIG. 1 shows a basic configuration of a Shack-Hartmann wavefront sensor 180. The 
Shack-Hartmann wavefront sensor 180 comprises a micro-optic lenslet array 182 that breaks an 
incoming beam into multiple focal spots 188 falling on an optical detector 184. Typically, the 
optical detector 184 comprises a pixel array, for example, a charge-coupled device (CCD) 
camera. The lenslets of the lenslet array 182 dissect an incoming wavefront and creates a pattern 
of focal spots on a charge-couple device (CCD) array. The lenslet array typically contains 
hundreds or thousands of lenslets, each on the size scale of a himdred nndcrons. Most Shack- 
Hartmann sensors are assembled such that the CCD array is in the focal plane of the lenslet array. 
[0007] A Shack-Hartmann wavefront sensor uses the fact that light travels in a straight line, to 
measure the wavefront of light. By sensing the positions of the focal spots 188, the propagation 
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vector of the sampled light can be calculated for each lenslet of the lenslet array 1 82. The 
wavejfront can be reconstructed from these vectors. 

[0008] To better understand one or more aspects of this invention, it is worthwhile to discuss 
the operation of a Shack-Hartmann wavefront sensor in more detail. 

[0009] In typical operation, a reference beam (e.g., a plane wave) is first imaged onto the array 
and the location of the focal spots is recorded. Then, the wavefront of interest is imaged onto the 
array, and a second set of focal spot locations is recorded. FIGs. 2A-F illustrate this process. 
[00010] Usually, some optical system delivers a wavefront onto the lenslet array which samples 
the wavefront over the tiny regions of each lenslet. The lenslets should be much smaller than the 
wavefront variation, that is, the wave-front should be isoplanatic over the sampled region. When 
the CCD array is in the focal plane of the lenslet array, each lenslet will create a focal spot on the 
CCD array. The location of these focal spots is the critical measurement, for this reveals liie 
average of the wavefront slopes across each region. That is, the shift in the focal spot is 
proportional to the average of the slopes of the wavefront over the region sampled by the lenslet. 
Software may compute the shift in each focal spot. 

[00011] If the wavefront is not isoplanatic, the quality of the focal spot erodes rapidly and it 
becomes more difficult to determine the peak location. 

[00012] However, where the isoplanatic condition is satisfied and where the focal spot shift is 
consistent with the small angle approximation of Fresnel, then the focal spot shift is exactly 
proportional to the average of the wavefront slope over the lenslet. 
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[00013] The incident wavefront is then reconstructed from the measurements of the average of 
the slopes for the hundreds or thousands of lenslets in the array. 

[00014] It was stated above tiiat the shift in the focal spot is proportional to the average of the 
slopes of the wavefront over the region sampled by the lenslet. This will now be shown as 
follows: 

[00015] We begin by examining the types of wavefronts incident upon a lenslet and calculating 
their spatial irradiance distribution in the focal plane. We then resolve what this distribution tells 
us about the average value of the slopes of the sampled wavefront. 

[00016] As noted above, most Shack-Hartmann sensors are assembled such that the CCD array 
is in the focal plane of the lenslet array. Therefore the appropriate complex amplitude 
distribution, Ff(u, v), of the field in the focal plane is given by the Fraunhofer diffraction pattern: 

[00017] 1) Fj{u, V) = -L- j jF^(x,>;)exp[-i|2(x« '^yv)y{x,y)dxdy 

—00 

[00018] where the pupil function describes the lens aperture. The function is 1 inside the lens 
aperture and 0 outside of the aperture. In the example that follows we will build the pupil 
fiinction using the rectangle function, which is defined as: 



[00019] 2) n(x)s< 



0 for \x\ > i 

1 for 

1 for \x\ < i 
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[00020] The field Fi(u, v)incident upon the lenslet is related to the wavefiront \|/ (x, y) incident 
upon the lenslet through: 

[000211 3) Fiix.y) = Gxp(ri\\fikx,ky)) 

[00022] where k is the wavenumber, Infk. The spatial irradiance distribution in the focal plane 
is related to the field amphtude by: 

[00023] 4) I(x,y) == Fix,y)F''ix,y) 

[00024] The first step was to model the portion of the wavefiront sampled by the lenslet. We 
begin with planar wavefironts and advance to wavefironts with more structure. Consider the first 
six simple cases in a Taylor series: 



[00025] 



[00026] '^Vioix.y) = X 



[00027] ^oi(x,y) =y 
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[00028] V|/2o(^,;^) = ^ 



[00029J "^w^^.y) = 



[00030] M/o2(^»>') = 



[00031] These six wavefronts were then analyzed using equation (2). Before presenting the 
results for the corresponding intensity distributions, a few functions must be defined. First is the 
sampling function 



sinx 

[00032] (5) sincx s 



[00033] then the error fimction: 



[00034] (6) erf{x) ^-jrje^ dt 



[00035] and the imaginary error fimction 



[00036] (7) erfi{x)s-ierf{ix) 
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[00037] and finally the exponential integral function 



00 



[00038] (8) eiix) = - j^^t 



[00039] The intensity distributions can now be written in terms of these general functions. The 
first is the well-known result: 

[00040] (9) Ioo(x,y) = (^'sinc^(^:c)sinc'(^y) 

[00041] The next two results are also for planar wavefronts wilh tilt. They are: 

[00042] (10) /^o(^>>') = {^""'^A^i^-^y^Afy) 



[00043] and 



[00044] (11) /oi(x 



[00045] We expect the symmetry under interchange of x and y in equations (10) and (1 1) since 
space is isotropic and the lenslet has a discrete symmetry under rotations by integer multiples of 
90°. 
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[00046] For the case of nonplanar wavefironts the results quickly become more complicated. To 
present these results, first define the intermediate variables 



[000471 (12) a = iZE^^t^ and P = ^2^^=^ 

2/A ^J'^ 



[00048] The intensity distribution can now be written as: 
[00049] (13) hi^i.^.y)'--]^ierfi{z^>^-^ 

[00050] For both l2o(x,y) and Io2(x,y) we see the symmetry in the solutions under interchange of 

X and y. The intermediate variables now become: 



1000511 (U) a = 22^ and P " 



[00052] and the irradiance distribution in the focal plane is: 



[00053] (15) hiix^y) = -i(|)\er/z(za) + erfiiz^yiierfiiizo^i^ er^^^^ 



[00054] Finally we consider the case where the input waveifront exhibits some torsion. First 
define the intermediate variables: 
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(000551 (16) a, = ?^!^ and «2 = 

(000561(17) P, = ^2^ ^ 152 = ^=^ 

[00057] Then the irradiance distribution can be written as: 
[00058] (18) I^^ix^y) = [-^^(eiiia^^^)-eii^^^ 

471 / 

X (ei(-/aiPi) --e/0'ai|l2) -eiO'ajPi) + efHajpj)) 

[00059] Clearly the solutions in equations (9)-(l 1) are even functions and they are separable. 
That is, the irradiance distributions for these three functions can be written X(x)Y(y). The next 
three fiinctions are not separable as the function arguments mix the variables x and y. However, a 
careful analysis of the constituent functions in equations (6)-(8) reveals that the combinations in 
equations (13), (15) and (18) are even functions also. 

[00060] Now we address the issue of how the amplitude of these input wavefronts affects the 
focal spot position. Consider that the wavefront has an amplitude k. The incident wavefronts 
become: 

[00061] ^ooi^yy) = K 
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[000621 ^io(x,y) = KX 



[00063] ^Qi(x,y) = Ky 



[000641 ^20^^>y^ ^ 



[000651 M'liUj') = W 



[000661 ^ 



[000671 These new wavefronts were again propagated through the lenslet using equation (1). 
As one might expect, there was no change in Ioo(x,y). The irradiance distributions for the odd 
parity terms become: 

[000681 (19) /io(^,>') = (^y^^^\^i^-^j)y^o\^y) 

[000691 and 

[000701 (20) loi y) = @'sinc'(^;c) smc\^(y - k/)) 
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[00071] For the higher order cases, tiiere is a slight rewrite. For example for l2o(x), the 
intermediate variables become: 



[000721 (21) a=?5^^ aad P = '"^^^ 



[00073] and the irradiance distribution in equation (1 3) is then multiplied by k"^'^ . For the 
torsion case y ii(x, y) the intermediate variables are: 



[00074] (22) ai 2A "2 2A 



100075] (23) .p,=2^H^^ and P2 = ^^^^^ 



[00076] and the irradiance distribution in equation 18 is then multiplied by k~^. 

[00077] However, the only peaks that shift position are Iio(x,y) and Ioi(x,y). And both of these 

shifts are exactly in the focal plane. The question now becomes how does this compare to the 

average value of the wavefront slopes? 

[00078] The average value of the wavefront slope is given by 

d/2 d/2 



[00079] (24) jn^ -4/2-d/2 — 



n 
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[00080] and similarly for my. The six input wavefironts and their average slopes m = (mx, my) 
are: 



[00081] M'oo(^'>') = K = (0, 0) 

[00082] ^ioix,y) =Kx = (0, 0) 

[00083] Woiix,y)=Ky m = (k, 0) 

[00084] ¥20 C^'^') = Kjc^ = (0,k) 

[00085] ¥11 C^,;') = Kxy = (0,0) 

[00086] Wo2iXyy) = = (0,0) 



[00087] So for these six cases, which are the most physically relevant, the average of the 
wavefront slope exactly describes the shift of the peak in the focal plane accounting for the 
distance offset of the focal length. Because of this exact relationship, it is convenient for us to 
define the focal spot location as being the location of the peak intensity of the focal spot. 
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[00088] For the purposes of this discussion, we define "isoplanatic" as the condition where the 
waveiBront is well approximated over an area the size of the lenslet, by a plane wave. 
[00089] This analysis buttresses the lore of Shack-Hartmann wavejfront sensing. The shift in the 
focal spot position is directly proportional to the average of the wavefront slope across the 
lenslet In the small angle limit inherent in the Fresnel approximation, this is an exact result for 
the isoplanatic cases and also true for the higher terms examined, although focal spot location is 
ambiguous in these cases as well as for additional higher order terms. An important issue that has 
not been discussed is how combinations of these terms behave. This is an important issue that is 
quite relevant physically. 

[00090] FIGs. 3A-F and 4A-F show some exemplary incident wavefronts and output 
irradiance distributions in the focal plane. The lenslet is a square of size d = 280 n and has a focal 
length f = 28 nmi; the wavelength X = 325 irni. 

[00091] FIGS. 3A-C show the functional form of some exemplary incident wavefronts which all 
have a generally planar structure. FIGs. 3D-F graphically shows the corresponding wavefronts 
incident upon the lenslet. Note that in these examples over 95% of the light is confmed to an 
area the size of a lenslet. The x and y axes are in microns. The vertical scale in the first colunm is 
also in microns. The vertical axis on the second column represents intensity and is in arbitrary 
units. 

[00092] FIGs. 4A-C show a series of increasingly more complicated incident wavefronts. 
Clearly the focal spots are degrading rapidly where the incident wavefront has a non-planar 
structure. Notice the concept of a peak is ambiguous here and a center of mass computation of 
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focal spot location could be highly problematic. Also, only a small portion of the light falls 
within the lenslet area. The x and y axes are in microns. The vertical scale in the first colunm is 
also in microns. The vertical axis on the second column represents intensity and is in arbitrary 
units. 

[000931 FIGs. 3A-F and 4A-F validate the well known mle of thumb: the wavefi:ont sensor 
works best in the isoplanatic limit. In other words, the lenslet size should be much smaller than 
the variations one is attempting to measure. Of course practical wavefront sensor design includes 
many other considerations too and one cannot shrink tixe lenslet size arbitrarily. But the message 
is that in the isoplanatic limit we expect to find well-defined focal spots that the center of mass 
can handle with higher precision. 

[00094] So, to recap, a Shack-Hartmann wavefiront sensor is a powerfiil tool for optical analysis, 
however the precision of the wavefiront data that it produces maybe less than is desired at times, 
and the speed and ease by which an incident wavefront may be reconstructed are also sometimes 
less than desirable. 

[00095] First, it is difficult to precisely locate the focal spots produced by the lenslet array. 
Although the center of mass may be computed as an approximation of the focal spot location, 
more precision is desired. The basic problem with center of mass algorithms lies with deciding 
which pixels to include in the calculation. All cameras have a background level of electronic 
noise. The center of mass algorithm is especially sensitive to noise away from the brightest pixel 
because of weighting of the calculation by distance. To exclude this effect, a threshold is usually 
applied that excludes all pixels below some absolute brightness level. This threshold may be 
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further refined after this preliminary step according to how bright the light is under a lenslet. To 
see the effect of this threshold, a calculation cau be done using the center of mass. The sinc- 
squared intensity pattern can be overlaid with a discrete section of CCD pixels. The continuous 
intensity pattern is broken iuto discrete pixelated values as counts. Then the center of mass 
algorithm can be applied to the pixelated values. Then a calculation can be done where the sinc- 
squared pattern is translated into small steps across the fixed CCD. This corresponds to the 
situation where tilt is added in small steps to the beam as it hits the lenslet. The center of mass 
location can be plotted against the known tilt added to the beam. When this process is done, it 
will be seen that the center of mass will follow a straight line up to the point where a new pixel 
starts to receive enough light to be counted (or another pixel got too little light.). At that point, 
there will be an abrupt jump in the line. The appearance will be a line with lots of little kinks in 
it. On the average, the input translation versus center of mass calculation will produce a line 
with a slope of one. So it is evident that if there are kinks in the line, the slopes of all the little 
straight sections of lines must be something different than the ideal slope of one. The kinks and 
the non-unity piecewise slopes indicates that center of mass calculations have a fundamental flaw 
that makes them unsuitable for being the basis of wavefiront measurements if extremely high 
accuracies are desired. 

[00096] Secondly, in general, a lenslet boundary will not line up neatly with a pixel boundary. 
Instead the boundary line will likely cross through pixels. For precise wavefiront measurements, 
it is important to precisely determine the location of the vertices or boimdaries of lenslets in a 
lenslet array with respect to pixels (e.g., CCD array) of a corresponding wavefront sensor. This 
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knowledge is needed in order to assign locations of the boundaries of the lenslets to locations on 
the CCD. 

[00097] Third, it is difficult to take a series of hundreds or thousands of measurements of the 
average of the slopes and reconstitute the incident wavefront. This process of taking the 
measurements and reconstituting the wavefront incident upon the lenslet array is called 
reconstruction and can be quite computationally intensive and time-consunndng. 
[00098] Accordingly, it would be desirable to provide a system and method of more accurately 
locating the focal spots in a wavefront sensor. It would also be desirable to provide a system and 
method for to precisely determine the location of the vertices or boundaries of lenslets in a lenslet 
array with respect to pixels (e.g., CCD array) of a corresponding wavefront sensor. Furthermore, 
it would be desirable to provide a system and method for reconstructing a wavefront from 
measurement data provided by a pixel array in a wavefront sensor. Even fiirthermore, it would 
be desirable to provide a system and method for efficiently computing polynomials for defining a 
wavefront detected by a wavefront sensor. And, it would be desirable to provide such a system 
and method which overcomes one or more disadvantages of the prior art. 
[00099] The present invention comprises a system and method for sensing a wavefront. 
[000100] In one aspect of the invention, a method of determining a location of a focal spot on a 
detector array improves upon a "pure" center of mass calculation. When a near plane wave 
illuminates a square lens, the focal spot has a sine-squared shape. This knowledge can be 
exploited to generate a nimiber of improvements to calculating the focal spot peak location. 
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[000101] Several techniques show improvement over pure center of mass methods for locating 
a focal spot. One technique uses the knowledge that the sinc-squared function has minimum in it 
before the side lobes appear. The focal spot location algorithm can be written so these minima 
are located, and then the pixels that are identified as being proper to include in the center of mass 
calculation will be those that lie within the minima. Anotiier technique is that the shape of the 
sinc-square function can be fit to values that show up in each pixel. This can be done with least 
squares method or by "sliding" methods. Then the equation for the sinc-squared gives the 
location of the peak. Still another refinement is to use a different weighting than a linear weight 
firom center in the calculation of center of mass. The effect is to reduce weight given to dimmer 
pixels. In particular, this weight can be matched to provide piecewise linear slopes of one, 
although this reduces the dynamic range of the sensor unless another method is found to 
eliminate the kinks. 

[000102] In another aspect of the invention, a method of determining a location and size of a 
lenslet with respect to a detector array exploits a-priori knowledge of the spacing of the lenslets 
in a lenslet array. Such a-priori knowledge can be produced by direct measurement of the lenslet 
arrays by microphotography, by knowledge of the photolithography process, or by other means. 
In a particularly advantageous manufacturing process, the spacing between lenslet is extremely 
regular. Accordingly, when one compares algorithms that calculate focal spot locations, the 
better algorithm is the one that comes up with a the locations in the most regular grid. This can 
be done on a theoretical basis using the sinc-squared patterns as inputs, or on experimental data 
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[000105] In still another aspect of the invention, a method of determining a set of polynomials 
defining a wavefront from a plurality of focal spots produced on a detector array, includes . 



BRffiF DESCRIPTION OF THE DRAWINGS 



[000106] FIG. 1 shows a basic configuration of a Shack-Hartmann wavefi*ont sensor. 
[000107] FIGs. 2A-F illustrate a reference beam and a wavefiront of interest being imaged 
onto a detector array. 

[000108] FIGs. 3A-F show some exemplary incident wavefironts and corresponding output 
irradiance distributions in the focal plane. 

[000109] FIGs. 4A-F show some additional exemplary incident wavefi-onts and 

corresponding output irradiance distributions in the focal plane. 

[000110] FIG, 5 shows simulated CCD data. 

[0001 11] FIG. 6 illustrates center or mass calculations. 

[000112] FIG. 7 plots error in the center of mass calculation. 

[0001 13] FIGs. S-24r 8-25 illustrate other aspects of the invention. 



DETAILED DESCRIPTION 



[0001 14] OPTIMIZED METHODS FOR FOCAL SPOT LOCATION USING CENTER OF 
MASS ALGORITHMS. 
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[000113] 




DETAILED DESCRIPTION 



[000114] OPTIMIZED METHODS FOR FOCAL SPOT LOCATION USING CENTER OF 
MASS ALGORITHMS. 

[000115] We want to know the location of the focal spot and we use the center of mass 
technique to approximate this location. The question is: how well does the center of mass 
approximate the focal spot location? To resolve this issue we begin with a special case, and then 
move to the situation encountered in the laboratory. 

[000116] Consider the special case of a spatial irradiance function g(x) being symmetric about 
the peak in some local region. In other words, the function has even parity about the peak xO in 
some local region which is smaller than the integration domain for the center of mass. The 
requirement for even parity can be written as: 



[000118] Since g(x) is even, the antiderivative of this function is odd on this same interval. 
Defining the antiderivative as: 



[000117] 



(25) 
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[000119] |g(x)dfi; - , (26) 



[000120] the parity condition becomes: 



[000121] G(x-Xo) - -G(sH-Xq). (27) 

[000122] Now we consider the center of mass computation. Defining the center of mass on 
some domain [a, b] as: 

fxg{x}ds 

[000123] {2S) 

[000124] we now have the most general formula. But because of the parity of g(x), we are 
motivated to consider a symmetric domain centered on the peak. This implies: 

[000125] XQ-ff = -(xQ-fi). (za) 



[000126] If we consider p to be radius we can define: 



[000127] P=xa-a (3§) 



20 



I 



Docket No. WFS.017 
PATENT 



[000128] and the center of mass in equation 28 now becomes: 



[000129] 




[000130] Since G(x) is odd over this domain, the integral vanishes. And using equation 27 the 
center of mass can be reduced to: 



[000132] So when g(x) is symmetric about the peak and the domain is centered on the peak, the 
center of mass is exactly the focal spot location. 

[000133] Due to the discreet nature of the CCD array, the pixel addresses are integers and this 
precludes one firom selecting integration boundaries that are centered exactly about the peak. This 
leads to the consideration of how the center of mass shifts away from the peak. The approach 
here is to consider the problem as having two parts: a piece with exact symmetry and the piece 
which shifts the center of mass away from the focal spot location. 

[000134] For the arbitrary domain [a, b], no longer centered about xq, the center, c, of the 
domain is: 



[000131] 
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[000135] c = (33) 



[000136] This leads to the definition of an asymmetry parameter t which measures how far the 
peak is from the domain center; 

[000137] T = jCq-c, (34) 

[000138] The next step is to defme the radius p of the symmetric piece: 

[000139] P = Mjn(jCQ - a, 6 - Xq) . (35) 

[000140] The center of mass can now be written in terms of the symmetric and asymmetric 
component as: 

[000141] 5 = ^g^P 



[000142] which reduces to: 
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[000143] 



X = 



(xo-p)G(xa + p) + (xQ + p-2T)G(xQ+p-2t)- j G{x)dx 



(37) 



[000144] As T -»0 tlie integration boundaries become symmetric about the peak and the center 
of mass approaches the pealc, i.e. x -^«o- 

[000145] The error 5 = x - xq in the center of mass result can be expressed in terms of the 
asymmetry parameter x which is a key to a subsequent correction scheme. For the case of a 
sinc2x function, this error is given as: 



[000146] 



S = 



k(p) - ]n(p - 2t) + ci(2(p - 2t)) - gi(2p) 



(38) 




,2 



,2, 



[000147] where ci(x) is the cosine integral fimction: 



[000148] 




[000149] and Si(x) is the sine integral function: 



[000150] 




D 
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[000151] We now present three correction strategies to improve the center of mass computation 
and recover a more accurate focal spot location. 
[000152] Integration Domain Boundary Method. 

[000153] One method for implementing the center of mass to exploit the results of the analysis 
above locates the brightest pixel under a lenslet and then marches out in the ±x and ±y directions 
to fmd the four minima. These four pixels define the boundary of a rectangular integration 
domain. The smallest of the four minima is used as a threshold and this value is subtracted from 
all pixels in the integration domain. 

[000154] FIG. 5 shows simulated CCD data for illustrating how an irradiance distribution is 
represented discretely. Here the incident wavefront was isoplanatic and with a small tilt 
component. The rectangular box defmed by the white lines represent the integration region. Here 
the error is pixels. 5*^ = (0.013, 0.002) 

[000155] FIGs. 6A-B shows that the center of mass calculation does not match the focal spot 
location in the general case (FIG. 6 A) due to an asymmetry. The center of mass can be thought 
of as having two components: a symmetric piece (light gray) and the asymmetric component 
(dark gray). FIG. 6B shows that removing the pedestal greatly reduces the shift caused by the 
center of mass computation by reducing the contribution of the asymmetric piece. The 
asymmetric piece is the contaminant that pulls the center of mass value away from the peak 
location. Successfid center of mass strategies will minimize the effect of this piece. 
[000156] FIGs. 6A-B illustrate the wisdom of the Integration Domain Boxmdary method. First 
of all, it endeavors to symmetrically bound the peak, which minimizes the error in the center of 
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mass computation. This also gobbles up more pixels than a thresholding method, which provides 
more information. Now the asymmetric contaminant is also minimized near the minimum. The 
total area in the darker region is what causes the error in the computation. Near the minimum, 
this piece has the smallest height. And after removing the pedestal (the dark count from the 
camera), the piece is even smaller. The width of the piece is now less than the width of a pixel. 
So this simple scheme exploits the lessons from the analysis above all without having to provide 
a threshold a priori. 

[000157] There is a drawback to this method in that as the distance between the asymmetric 
piece and the peak increases, the lever arm also increases. In other words, a unit area added to the 
center of mass is more harmfiil at larger distances. However it appears tihiat in the cases of square 
lenslets with well-defmed suic2x sinc2y focal spots, the lever arm effect is more than 
compensated for by the height minimization at the minimum. This seems to suggest against using 
the second minimum as the lever arm doubles. 

[000158] Some basic numerical studies were on simulated camera files created using the 
propagation equations defined above. These data did not include either noise or cross-talk from 
neighboring lenslets. The advantage of simulated data is that one has exact knowledge of the 
peak and can exactly quantify the error in the center of mass calculation. For comparison 
purposes, a classic center of mass with a threshold of 12% was used. The average error 
5^ = (5j^, sp for each method was then computed for 4096 samples all where the peaks were 
within one pixel of the projection of the lenslet center. The standard deviation of the errors are: 
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[000159] 0/ = {0.015, 0.015 } pixels for the pure center of mass method, and 



[000160] <yfi^ = {0.010, 0.01 1 } pixels for the integration domain boundary method. 
[0001 61] Modified Integration Domain Boundary Method. 

[000162] The above Integration Domain Boundary method is an excellent starting point for 
improving the center of mass calculation. In the cases where the form of the irradiance 
distribution is loiown, an additional improvement is possible. Since the shape of the distribution 
is known, one can look at the difference between the center of mass and the asymmetry 
parameter and make a quick correction based on a two-dimensional version of equation (38). In 
the sample environment tested, the additional correction was able to remove 90% of the error, 
reducing the error by an order of magnitude. 

[000163] FIG. 7 shows the error in the center of mass computation as a function of peak offset 
from the center of the pixel. The axes describe the distance the peak was offset from the center of 
the pixel. The color of the block reveals the magnitude of the error. As expected, there is no 
error when the peak is centered in the pixel. The maximum error in this example was 0.03 pixels. 
The effect of the discrete change in integration boundaries manifests distinctly as an abmpt 
change in the error. 

[000164] The error variations in FIG. 7 were studied in the x and y components. The next step 
was to plot the error components as a function of the asymmetry parameter. A surprisingly 
simple variation was observed and is shown in FIG. 8. FIG. 8 plots the center of mass error as a 
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function of the asymmetry parameter. Here the graph shows the x component of the error plotted 
in FIG, 7. Note that the functional relation between error and the asymmetry parameter is linear 
with two distinct branches. The errors are essentially linear with respect to the asymmetry 
parameter x with two distinct branches. This implies that a simple linear correction can be used to 
push the center of mass closer to the focal spot location. For these data the results of the linear 
fits were: 

[000165] positive branch error = (0.043756 ± 0.000004) x + (0.00000 =b 0.00002) 
[000166] negative branch error = (-0,0305515 ± 0.000004) x + (0.00000 ± 0.00002) 
[0001 67] The beauty of this method is that the slopes for the positive and negative branches are 
functions of the wavefiront sensor: As such, they only need be determined one time. The method 
is to compute a center of mass and an asymmetry parameter. Then determine which branch of the 
correction is appropriate, and apply a simple linear boost to nudge the center of mass closer to the 
peak location. The nudge is apphed separately to the x and y components of the center of mass 
prediction. 

[000168] FIG. 8 allows one to see graphically how precise this correction can be. The vertical 
spread in each cluster of points is approximately 0.002 pixels. The linear correction basically cuts 
this in half, leaving us with an error of approximately 0.001 pixels. Attempts to reduce the error 
further by handling the problem in two dimensions failed. 

[000169] For example, if the asymmetry parameter was at the maximum oft = 0.3 pixels, the 
correction is: 

[000170] (0.043756 ± 0.000004)0.3 + (0.00000 ± 0.00002 ) = -0.013 13 ± 0.00002 pixels. 
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which is added to the computsd center of mass. This will bring fee computatLon to within pixels of flie true 
focal spot locaticdL 

43 FEiiiinI pixels 

AnoliiEa' strataegy has been davebpad in callabciiatian with J. Rdler. This is an itea'ativd Bchjeme. The gist cf 1h& 
metod is that ana computes the center of ihass and d&velops an apprQximate' focal spot location^ (%q, ) . This looa- 
tion shmalfl ftgi yi ly h p^ mtiun a pixel cf fee trjg fjcal locaH u n, A xic w liiLc gration danuin is now defined, centered 
on thfl centJar of massy alua The size cf the doraain is given by the Tsveticov criteria: use all ftis data within the first 
mniina. We know ftoni equation 9 1hat the half^ p = ^/rf. So tiaemtegration domain is [x^-p^x^+p] and 

the inflegration rang^ is [y^-p^yQ + p] - Using these new boundarissj a new cento of mass is can^juted, (jc j, * 
along witli new boundaries. The picxsess continuBs until the distance betweenpeak locations on consecative itEsralions 
is smalleir than same arbitrary acceptance pararaeto^ 

Ttie problem is fliEft now the integration boundaries are no bnger integer pixel values, but they lie within the pixels. 
So tiie heart of this schjeme is to find a way to ^portion the ligiit; within a pixel and maintain the spe^ advantaga of 
ttie cEnter of mass. While Ihe encact counputafciGin of ftus apportionment can he messy, the exact solution can be 
Epprosimated quite wdl with a polynamkL 

Oonaider the cases of sinc^ in one dimension. The parameter 9 describes vtoere the rnicdmum fells within the 
pixel Par example, if the ininimum is in the middle of the phcdj <p - 0.5 ; if the minimum is three quartss of the 
way dsmm the plxd then <p «- 0.75 . The function is a weighting functioai varying from zero to one which 
describes how much of the light to include from the pixel. In this case the weighting function is given exacfly as 

[000171] 
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A 



(41) 



Itie constiuit arepieassats 



a » JS^(2ic)<^ 141SlSlS7613.2fi2S45Q24S7801S229g7494!29142433B492dS0..« . 



(42) 



Notice ffiat this is far pix&ls on the riBht-hflnd side. One. ssm exploit ftie evea paiily of fc© anicAc fimctioEi to get Hie 
valUiBS for the left-hand Bide. 

The pmtfleia is:of course that a naive implemoatattei ci Si^} m][ slow down the algorithm ttCT:iqndau5ly. W can 
es^jloit ihe fact th^jt Qver to doniam iatearest; ^ € [Oj,!]^ te funcfion^ a rather Bimple balTairiG(c FigiiiB 7 heJow 
shoTO ths A»si^1fcg fimctdmi as a solid curve. The poirils arethe results aEpiolyrwfijQial least ^iiares-fittkroiagji (st^ 
BiDfii. THuBi agEeeittBilt is qiaitegooi 

We caiL write &e kppicxdznate solution fl(9) as 



(43) 



where n lepresmls ftie dggree of fit In figure 7, n - 9. Nctice thatthis fbrmulatiQn avdds the amhiguoiLs qnaiMly (P. 

Of CQiaiaa. a diiectimplBittetita&m of Bqiiafon43 wduM be ponderOuily dow^^^ erode the speed advantage of Ihe 
center ofraass tediidquiB. We shew tTO cH^^ enjcodeftiEBBiesidtsinOH-: 



dQUbla c [lo] f 



c[o] 
cm 

c.[63 
[000173] 



= 7a3-73205Sa9474:41f 
= -3448.7204133862506; 
=: S172 .400852 00703 ; . 
= -573g-32«837D7a93gj 

= S(S45,aflaa7a74&s38fl| 

« - ©33 . i.6021^f S4S233ia j 
« 13a , BQaSf 2.97iSf57St52| 



//DBF fit GbQfficients 
//constant terin 
//linQscr term 
//coaf floiiarit f or 



[000174] 



Figure 7: The wsightrng fbncticin <B^{<X}) for partial pixels. If the mTnimuni falls v^-fliia the pdxcl at some point q), then the 
porticm nf h^t to include in the ocaatEx of maas computatim ia fii[(<p) ^ The cmvc is the exact salutkm in cquatiaa41 aadtiM 
doia aic the naulta of a ninlii'Orier polyninTnial Icjat aqugnea fit 
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C[7] = -9 •792719111240226 J 
C[B] = 0.29071£557781€4:25; 
C[9] = -0.001729454319027303; 

doublQ OKiaga {double * double phi, Int n) 

// c iiSF fit cosffiGientB for eguatian 41 
// phi pixel fraction 

// n degree of polynomial least squares fit 
// iBf weighting factor to use Iomega [phi] ) 
{ 

double Isf = Qji 
for (int 1 = n; i > Of i--) 
Isf = phi * (Isf + c[i] ) f 

iBf 4= C[Q] ; 

return Isf j 

} 



//coefficient for x*9 
//weighting function 



//Lap prediction 
//loop through oorcierB 
//build up answer 

//add In constant term 

//return the ensver 



This is an elegaatim] 
Bowerm^ a somewhat fastar msthod is 



n. 



double f (double ♦ c, double phi) //rapid polynoinial 
{ 

double Isf = c[0] -i- phi * (ctH 4- phi * {c[2] + phi * (c[3] 

4 phi * (Gt4] 4 phi * (ct5] 4 phi * (ct6] 

4 phi * (c[7] 4 phi * {G [8] 4 phi * (c [9] ))))))))) 



return Isf i 



The first metibod took 1 6% bngar to execate than the second molhod* A straightforward impletmentakLoa of equation 
43 took 13i times long^ than the seccDid method. 

S. SUMMARY 

We have shown farnmlly the w&ll-known result that the sMft in. ttio focal spoft location is propartional to tiie aveiago 
wflveftont stopa aDKSs a lensleL The result is exact in the isoplanatic and smallangb Hraits* Wsthm showed therela- 
tkmship betwean the fbral gpot location and the Gentair of mass caraputatioa The accuracy of the assumption that the 
centEsr of mass is tie fbcal spot locaitfon depends dii^ctiy vpm honsr symmetrically flTB mtegraticm domain heunds the 
peak. When the domain is axactly syinmBtric^ the center of mass exadly represents tfaepeak Of couise wittirteai dala^ 
the intsajgratian Ixiundanies are not centered upon 1he> peaks. dosed with three strate^gies Is impoiove the accuracy of 
the center of mass computatian and preserve the tone advantage of the meithod. 
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[000175] PRECISION MEASUREMENTS OF A LENSLET ARRAY. 
[000176] As discussed above, a critical step in data processing involves fixing the focal spot 
location and a common method used for quickness is the center of mass. But naive applications 
of this method have a precision of 0.01 pixels, which forms an inherent precision boundary for 
all computations. 

[000177] While the error in the center of mass technique is highly systematic, its manifestation 
in an experimental measurement is typically, to first order stochastic. The challenge then 
becomes how to remove these stochastic errors. The first tool in such situation is the least 
squares fit which is a simple way to remove stochastic errors. 

[000178] Here we describe how a least squares fit is used to precisely measure the size and 
location of a lens-let array. Here we assume a rectangular lenslet and sUde the grid around to find 
the best location assuming the focal spot locations are under the center of each lenslet. 
[000179] The least squares fit. 

[000180] The first decision is which parameters to extract fiom the least squares fit. Although 
we could assume a square lenslet configuration, we used a rectangular configuration. Of course 
we expect to find a square, but by solving for the two widths (dx, dy) independently, we are able 
to check the quality of the residt by seeing if the aspect ratio is unity. Also, this formulation 
allows for the more general case of rectangular lenslets. The other free parameter is the location 
of the lenslet origin, O^ = (xo, yo). 
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2.1 The problem 

The fitting strategy is elementary. Consider two rows of points, the points corresponding to the focal 
spot locations under each lenslet. We would want to find the height of a line which best bisects the 
two rows. 



(1) 



Here r is the row index and c is the column index and y is the abscissa of the focal spot location. The 
summation is across all columns and involves and arbitrary row r and the next row, r + 1. Clearly, if 
all the points rested upon two horizontal lines, the merit fimction % would be zero when h is the 
bisector. 

This method is easily extended to the other dimension by finding the line which best bisects each 
column. To formulate the full problem in two dimensions, we need to define a few quantities. We 
already have the origin vector O, 



O = 



and the lenslet width vector d, 



We need to define an position vector £, which specifies the row and column the lenslet: 



(2) 



(3) 



(4) 



Finally we note that the focal spot location for the lenslet in column c and row r of the lenslet array is 



(5) 



so we are loading the x and y locations into a vector for efficient computation. We note that the k vec- 
tor describes the number of rows and columns in the lenslet array that are projected onto the CCD 
array: 



k = 



(6) 



Finally we need a collection of basis vectors. In two dimensions, it is a spinor (Weisstein, p 1703) 



a = 



(7) 
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This structure points to the neighboring object. Although this is a two-dimensional problem, this 
methodology can be trivially extended to an arbitrary number of dimensions and it reduces the com- 
plexity of computer code. 
If we use // as a spacetime index to sum over the dimensions, the full merit function is 

^'=Z S (^.^VVk^/p.a/. (8) 

This problem is solved in the canonical fashion. All first derivatives of the fitting parameters are 
forced to zero simultaneously by solving a set of linear equations. That is, the foxir simultaneous 
requirements that 



9oX =0, 



(9) 



and 



=0, (10) 
generate a series of matrix equations. The solution to these equations follows. 
2.2 The solution 

First, we define some intermediate variables. There are two summation variables, n and a , They 
are given by 



1.1 



(11) 



and 



< = E^- (12) 

1. 1 

If every lenslet registers a focal spot, then these quantities are (CRC Standard Mathematical Tables, 
p54) 



a = 



2 



(13) 



(14) 



and 
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i 



2 

a = 



MV^)(2^r+l) 



(15) 



Finally, we define a determinant, A ; 



A 2 / I 2 



(16) 



We can now begin writing out the solutions. The lenslet origin is 



O = — P°li. 11 



2A.. 



(17) 



The lenslet size is 



2A.. 



(18) 



We now t\irn our attention to the computation of statistical uncertainties. The first quantity to com 
pute is the sample standard deviation (Bevington, p 19). 



(19) 



This allows us to compute the uncertainties in the size of the lenslets and the position of the array. 
The uncertainty in the lenslet origin is 

(20) 



and the uncertainty in the lenslet width is 



(21) 



3. Example 

A sample lenslet and camera were chosen from the Columbus experimental configuration. This sen- 
sor uses a 60 X 60 array of square lenslets. Each lenslet is 280 across with a spherical profile and a 
focal length of 28 mm. The array is etched into fused silica and the lenslets have a 100% fill factor. 
The camera used was an SMD model 1M15, with a 14 ^ pixels in a 1024 X 1024 grid. The output is 
digitized to 12 bits. The lenslet array was placed approximately 36 mm (1.3 focal lengths) from the 
CCD array. 



EiMibraiyMensletMenslet OS.fin 



34 



17^)5:23, enm 



I 



3.1 Single measurement 

By design, the lenslet array is much larger that the CCD array, so a subset of the lenslets are used. 
Here we used a 50 X 50 subset of the lenslets. For these 2500 lenslets, the array, was found to be 

size = 20.0005 ± 0.0008 x 20.0025 ± 0.0008 pixels, 
and the location of the origin of the array was found at the pixel address of 

location = (16.92 db 0.02, 15.58 ± 0.02) pixels. 
The aspect ratio of the lenslet is then 

"^=i^?fS 

implying that the lenslets are very close to being square. 

Of course this measurement does not simply measure the regularity of the lenslet array, it also 
measures the regularity of the CCD array. But the basic technique demonstrates clearly that an 
instrument with an inherent low precision limit can be an effective tool to make precise measure- 
ments. 

3.2 The tilt experiment 

The data above are part of a larger set of data collected by T.D. Raymond. An optical fiber was set 
750 mm from a coUimating lens and a shear plate was used to check collimation. These data are the 
baseline file, xOyO.bvf. Then a differential micrometer was used to translate the fiber in a plane per- 
pendicular to the direction of propagation. This method introduced up to 400 ^r of tip and tilt in 50 ^r 
increments. The data were analyzed for each position and a subset of the results are shown in table 1 
below, 

A close look at the data reveals a very shght asymmetry, at the fringe of the experimental uncer- 
tainty. Most formally this number is ratio of the aspect ratio of the lenslet to the aspect ratio of the 
CCD array. The natural question becomes is the asymmetry in the lenslet array or the CCD array? 
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File origin 

xOyO 16.92 ± 0.02, 15.58 ± 0.02 

xOySO 16.79 ± 0.02, 15.58 ± 0.02 

xOylOO 16.67 db 0.02, 15,58 ± 0.02 

x0yl50 16.55 ± 0.02, 15.58 ± 0.02 

x0y200 16.41 ± 0.02, 15.58 ± 0.02 

x0y250 16.29 ± 0.02, 15.58 ± 0.02 

xOy300 16.16 ± 0.02, 15.58 ± 0.02 

x50y0 16.92 ± 0.02, 15.70 db 0.02 

x50y50 16.80 ± 0,02, 15,70 ± 0.02 

x50yl00 16.67 ± 0.02, 15.70 ± 0.02 

x50yl50 16.55 ± 0.02, 15.70 ± 0.02 

x50y200 16.41 ± 0.02, 15.70 ± 0.02 

x50y250 16.29 ± 0.02, 15.70 ± 0.02 

x50y300 16.17 db 0.02, 15.70 ± 0.02 

xlOOyO 16.92 ± 0.02, 15.82 ± 0.02 

xlOOySO 16.80 ± 0.02, 15.82 ± 0,02 

xlOOylOO 16.67 ± 0.02, 15.82 db 0.02 

xl00yl50 16.55 db 0.02, 15,82 ± 0.02 

xl00y200 16.42 ± 0.02, 15,82 ± 0.02 

xl00y250 16.29 ± 0.02, 15.82 ± 0.02 

xlO0y300 16.16 ± 0.02, 15.82 ± 0.02 

xl50y0 16.92 ± 0.02, 15.95 ± 0.02 

xl50y50 16.79 ± 0,02, 15.96 db 0.02 

xl50yl00 16.67 ± 0.02, 15.95 ± 0.02 

xl50yl50 16.55 db 0.02, 15.96 ± 0.02 

xl50y200 16.42 ± 0.02, 15.96 db 0.02 

xl50y25b 16.29 ± 0.02, 15.96 db 0.02 

xl50y300 16.18 ± 0.02, 15.95 ± 0.02 

x200y0 16,92 ± 0.02, 16.10 ± 0.02 

x200y50 16.80 ± 0.02, 16.10 ± 0.02 

x200yl00 16.68 ± 0.02, 16.10 ± 0.02 

x200yl50 16.56 ± 0.02, 16.10 ± 0.02 

x200y200 16.42 ± 0.02. 16.10 db 0.02 

x200y250 16.29 ± 0.02, 16.10 ± 0.02 

x200y300 16.18 ± 0.02, 16,10 ± 0.02 

x250y0 16.93 ± 0.02, 16.23 ± 0.02 

x250y50 16.80 ± 0.02, 16.23 db 0.02 

x250yl00 16.67 ± 0.02, 16,22 ± 0.02 

x250yl50 16.54 ± 0.02, 16.22 i 0.02 

x250y200 16.41 ± 0.02, 16.23 ± 0.02 

x250y250 16.30 ± 0.02, 16.23 ± 0.02 

x250y300 16.17 ± 0.02, 16.24 ± 0.02 

x300y0 16.93 ± 0.02, 16.35 i 0.02 

x300y50 16.80 db 0.02, 16.35 ± 0.02 

x300yl00 16.69 ± 0.02, 16.35 ± 0,02 

x300yl50 16.56 ± 0.02, 16.35 db 0,02 

x300y200 16.42 ± 0.02, 16.35 ± 0.02 

x300y250 16.30 ± 0.02, 16,35 ± 0.02 

x300y300 16.17 db 0.02, 16.35 ± 0,02 
Table 1 



size 

20.0005 ± 0.0008 x 20.0025 db 0.0008 

20.0006 db 0.0008 x 20.0026 ± 0.0008 
20.0005 ± 0.0008 x 20.0027 db 0.0008 

20.0005 ± 0.0008 x 20.0026 ± 0.0008 

20.0007 ± 0.0008 x 20.0026 db 0.0008 

20.0006 ± 0.0008 x 20.0027 i 0.0008 

20.0006 db 0.0008 x 20.0027 ± 0.0008 

20.0007 db 0.0008 x 20.0030 ± 0.0008 
20.0005 ± 0.0008 x 20.0029 ± 0.0008 

20.0005 ± 0.0008 x 20.0029 ± 0.0008 

20.0006 ± 0.0008 x 20.0029 ± 0.0008 

20.0008 ± 0.0008 x 20.0030 ± 0.0008 

20.0007 ± 0.0008 x 20.0029 ± 0.0008 

20.0006 ±0,0008x20.0028 ±0.0008 ' 

20.0007 ± 0.0008 x 20.0032 ± 0.0008 

20.0005 ± 0.0008 x 20.0030 ± 0.0008 

20.0006 ± 0.0008 x 20.0032 ± 0.0008 
20.0006 ± 0.0008 x 20.0030 ± 0.0008 

20.0006 ± 0.0008 x 20.0030 ± 0.0008 

20.0008 ± 0.0008 x 20.0032 ± 0.0008 

20.0007 ± 0.0008 x 20.0033 ± 0.0008 
20.0006 ± 0.0009 x 20.0028 ± 0.0009 
20.0006 ± 0.0009 x 20.0029 ± 0.0009 
20.0006 ± 0.0008 x 20.0031 ± 0.0008 
20.0006 ± 0.0008 x 20.0028 ± 0.0008 

20.0006 ± 0.0008 x 20.0027 ± 0,0008 

20.0007 ± 0.0008 x 20.0028 ± 0.0008 

20.0004 ± 0.0008 x 20.0028 ± 0.0008 
20.0007 ± 0.0009 x 20.0022 ± 0,0009 

20.0006 ± 0.0008 x 20.0022 ± 0.0008 

20.0005 ± 0.0008 x 20.0022 ± 0.0008 
20.0005 ± 0.0008 x 20.0022 ± 0.0008 

20.0007 ± 0.0008 x 20.0022 ± 0.0008 
20.0007 ± 0.0008 x 20.0022 ± 0.0008 

20.0005 ± 0.0008 x 20.0020 ± 0.0008 
20.0007 ± 0.0008 x 20.0018 ± 0.0008 

20.0006 ± 0.0008 x 20.0017 ± 0.0008 

20.0007 ± 0.0008 x 20.0019 ± 0.0008 

20.0007 ± 0.0008 x 20.0020 ± 0,0008 

20.0008 ± 0.0008 x 20.0019 ± 0.0008 
20.0006 ± 0,0008 x 20.0018 ± 0.0008 

20.0006 ± 0.0008 x 20.0018 ± 0.0008 

20.0007 ± 0.0008 x 20.0021 ± 0.0008 

20.0006 ± 0.0008 x 20.0021 ± 0.0008 

20.0004 ± 0.0008 x 20.0020 ± 0.0008 

20.0005 ± 0.0008 x 20.0020 ± 0.0008 
20.0005 ± 0.0007 x 20.0022 ± 0.0008 

20.0007 ± 0.0008 x 20.0021 ± 0.0008 



aspect 
1.00010 ± 0.00006 

1.00010 ±0.00006 

1.00011 ± 0.00006 
1.00011 ±0.00006 
1.00010 ±0.00006 
1.00010 ± 0.00006 

1.00010 ± 0.00006 

1.00011 ±0.00006 

1.00012 ± 0.00006 
1.00012 ± 0,00006 
1,00011 ±0.00006 
1,00011 ± 0.00006 
1.00011 ± 0.00006 

1.00011 ± 0.00006 

1.00012 ± 0.00006 

1.00013 ± 0.00006 
1.00013 ±0.00006 
1.00012 ±0.00006 
1.00012 ± 0.00006 

1.00012 ±0.00006 

1.00013 ±0.00006 

1.00011 ±0.00006 

1.00012 ± 0.00006 
1.00012 ±0.00006 
1.00011 ± 0.00006 

1.00011 ±0.00006 
1.00010 ±0.00006 

1.00012 ± 0.00006 

1.00007 ± 0.00006 

1.00008 ± 0.00006 

1.00009 ±0.00006 
1.00008 ±0.00006 
1.00007 ±0.00006 
1.00007 ±0.00006 
1.00007 ± 0.00006 
1.00006 ± 0.00006 

1.00005 ± 0.00006 

1.00006 ±0.00005 
1.00006 ± 0.00005 

1.00005 ±0.00005 

1.00006 ±0.00005 

1.00006 ±0.00005 

1.00007 ±0.00006 

1.00008 ± 0.00006 
1,00008 ±0.00006 
1.00007 ± 0,00006 
1.00007 ±0.00005 

1.00007 ± 0.00005 

1.00008 ±0.00005 



20.0007 ± 0.0008 x 20.0022 ± 0.0008 
The results of many measurements. A differential micrometer was used to introduce tip 
and tilt in 50 increments. The computations in equations 17-21 were then performed on the 
new data set. 
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Docket No. WFS.017 
- PATENT - 

[000181] WAVEFRONT RECONSTRUCTION FOR THE SHACK-HARTMANN 
WAVEFRONT SENSOR. 

[000182] As discussed above, one uses the results of myriad averaged slope measurements to 
reconstruct the wavefront incident upon the lenslet array. It is common to hear people describe 
reconstruction as an integration process. This is malapropos and ironic. The irony is that Nature 
has aheady done the integration for us; as shown below the measured quantity is already an 
integrated form of the incident wavefront. 

[0001 83] Given a wavefront \|/ (xl , x2) incident upon a lenslet, what exactly does the sensor 

measure in the isoplanatic and small angle limits? It measures the average of the wavefront 

slopes over each lenslet. So for an integration domjun of xi e [a, b], X2 e [c, d], tibie average 

of the \x derivative of \|/(xl, x2) is given by: 

[000184] 

[000185] 
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Figure/: The Shack-Hartmann sensor tracks a series of focal spots. This idealization shows the proems of measuring a spher- 
ical^ve on a sensor with just 16 lenslets. The first step, as represented by the first column, is to meWre a plane wave and 
ijjeasure a series of focal spot locations which are used as a reference. The next step, as depicted in the s^pnd column, is to 
introduce the wavefront of interest and determine the focal spot shifts. The critical information at this stageVfhe analysis is 
the focal spot shifts for every lenslet. 
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where the operator 3^ denotes partial differentiation along a line parallel to the \i axis. Defining the antiderivative of 
as 

the average of the slope is written as 

= A''\^^ia,c)-^^ia,d)-^^{b.c) + ^f^ib,d)) (3) 

where the domain area has be recast asA=-{d- c){b -a) and v represents the orthogonal direction. The message 
of equation 3 is quite clear: the measurement of the focal spot shift tells us about differences in the antiderivatives 
evaluated at the lenslet vertices. 

At this juncture, the reconstruction will typically follow one of two paths: either modal or zonal reconstruction. An 
example of the two different types of outputs is shown below in figure 2. Zonal methods output a value for \|r at each 
zone. In the Shack-Hartmann case, the wavefiront values are computed at the vertices. In the Southwell case, the 
wayefront value is computed for the center of the lenslet. The modal reconstructor output the amplitudes of some 
arbitrary set of functions used to describe the wavefront. A practical set of functions are (he Taylor monomials. With 
these amplitudes, it is trivial to compute any other rational function set. Of course, one can not easily convert the Tay- 
lor amplitudes into amplitudes for a irrational function set, such as one contaming Gaussians. 

Consider a fit to order d using a sensor with an ^ x Tj lenslet array. The output from the zonal reconstruction will be 

+ if (11 + wavefront values for each of the vertices, reducing the wavefront to a discrete data set. The modal 
reconstructor will yield ( J + 1 )(flf + 2)/2 amplitudes, and the wavefront is now in a continuous representation. As an 
example, consider a 6th order fit with a 50 x 50 lenslet array. The zonal fit will yield 2601 numbers; the modal fit will 
produce 28 numbers. So the information content of the zonal reconstraction is typically richer. One may expect to pay 
a computation burden for this extra information, but this is not necessarily the case as we shall see. 

Figure 2 below shows the reconstruction process schematically, A wavefront is imaged onto the lenslet array creat- 
ing a pattern of focal spots on the CCD array. The software then computes how far the focal spots have moved from 
the reference positions and creates a series of offset measurements, 5. The focal spot shift is related to the lenslet focal 
length /and the average wavefront slope m by 

5^/7= (4) 
Vl 

when the CCD is in the focal plane of the lenslet array. In practice, the radical in the denominator is so close to unity 
that it is ignored leading to 

8 = > . (5) 

The solution to the full form is 




(6) 



and the approximation is 



« = f (7) 
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643=0 



C4i = 0 
C4o = 0 



'zonal output 



\ 

modal outputs 



Figure^ Reconstruction inputs and outputs. In this example, an eigenstate wavefront of the Zemike polynomials, 2^, has 
beea/fmaged onto a 48 X 48 lenslet array. The focal spot shifts are computed and passed to either the zonal or the mb 
Dnstructor. A san^le output from each reconstmctor is shown. The zonal reconstructor provides a value for the wavefr^ 
t every lenslet vertex; the modal reconstructor provides an amplitude for every function in the basis set. 
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For example, an array with a focal length of 20 mm is imaged onto a CCD array with 24 ^im pixels. A shift of two 
pixels is measured. The wavefiront slope then is 1 .2 mr using the approximation in equation 7, and is 1 .200000864 mr 
using the exact form in equation 6. So for practical purposes, the approximation is adequate for lenslet arrays with 
long focal lengths. 

The measurement of the focal spot shift 8 then is a measurement of the average slope over the lenslet. In the 
isoplanatic limit the average slope can be thought of as a tilt. The effect then of a measurement is to reduce the wave- 
front into a series of planar patches of the size of each lenslet, with a tilt in each direction corresponding to the focal 
spot shift in that direction. 

To delve deeper into the details of the calculation, we need to introduce some notation. There are a series of iST off- 
set measurements where 7/ = is the number of lenslets. Eventually, the linear algebra of zonal reconstruction is 
going to force the lenslets to be ordered in a linear fashion. In the interest of expediency, we will use the same order- 
ing for both the zonal and the modal cases. 

The ordering of the lenslets is completely arbitrary: it can be row major, column major or even random. The posi- 
tion vector describes the linear ordering of the addresses. A sample p vector is 



p = ((1,1), (2,1), (3,1)...) 



(8) 



where the address pairs (u, v) denote rows and columns. On an ^ xt) array, the address (u, v) would be at position 
m(^-- 1) + y m the /7 vector. The reverse look-up procedure is a bit more involved. It is simple to compute, but in 
practice, it is convenient to construct a lookup table q which is really an ^ x matrix. A sample q matrix is shown 
below 



1 

r| + l 



2 

j\ + 2 



n 

... 2ti 



.(^-l)Tl + l(^-l)Tl + 2 ... ^Ti, 



(9) 



3. MODAL RECONSTRUCTION 
Equation 3 has already paved the road for modal reconstruction. Before beginning with the details of reconstruction, 
it is helpftil to defme the box operator □ for an arbitrary function over the domain jc^ e [a, 6] , g [c, d] . This 
operator combines the values of a function at the domain comers in the following way; 



□VW = V(a,c)-\|r(a,flf)-Y(6,c) + \|f(b,<i). 



(10) 



Having seen the box operator in an explicit form, we will now generaUze its formulation. Consider the rectangular 
domain show in figure 3. The vertices are labelled by a, P, y,and 8 where each vertex has an address (xi, xj). These 
vertices can be arranged in a vertex vector K written as 



V = 



(U) 



When a function operates on this vector, it returns a vector of functions. Explicitly, 



¥(a) 

V(P) 
¥(Y) 
L¥(8). 



(12) 



(To save space, subsequent vectors will be written in the transpose form.) This leads to a generalization of operators 
which act upon the vertices. Think of an interaction matrix T which describes the scalar interaction of the vertices. 
For the box operator 
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Figure 3: The vertices of a lenslet. Here <h^^^'e^tex points are labelled with the Greek letters a through 5. Each of these letters 
is really an address pair. The indexing of the vertices is tied to the indexing for the lenslet. If this is lenslet (u, v), then a is ver- 
tex (u-l, v-i;^P is (u, v-V, Yis (u, v), and Sis (u-l,v). 



-10 0 0 
0 10 0 
0 0-10 
0 0 0 1 



(13) 



and we can now write 



With this box operator, the average of the \i slope over the lenslet at address is expressed as 



(14) 



(15) 



In this case p is the predicted value. But the actual measurement is of 5^ . To find the wavefront which best 
describes the mea^emfents, one simply needs to minimize the merit function ' 

li = 1 



2 



where a,- represents an error estimate. However, since these errors do not depend upon the amplitudes, they can be 
treated as weighting factors and are withheld from explicit evaluation. 

All that remains now is to chose a set of basis functions. This set should be complete. Orthogonal sets are attrac- 
tive, but often involve imneeded overhead as shown in [3]. The set we will choose is the basis set for all rational func- 
tions: the Taylor monomials. These monomials are given by 



and the wavefront described as a tfth order expansion would be written as 

d I' 



(17) 



(18) 



We specify this set because it contains all the information needed to construct the functions of jail other rational basis 
sets, such as the polynomials of Zemike or Legendre. It is easiest to work in the Taylor basis and then use an affine 
transformation on the Taylor amplitudes to compute the amplitudes in the other bases. These transformations matri- 
ces are written in terms of integers, so there is no loss of precision. 

Also, the Taylor basis far easier to work with. Consider the case of the antiderivatives. They can be written to arbi- 
trary order as 
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(19) 



and 



^2 \^V^2> T+T I 2 



(20) 



By comparison, writing out the two antiderivatives for each of the 66 Zemike polynomials through 10th order would 
fill several pages. 

Anyway, the following analysis does not depend upon the basis set. The minimization of the merit function pro- 
ceeds in the canonical way: the derivatives with respect to each amplitude are set to zero. 



(21) 



This will generate a series of « = + 2)/2 linear equations to solve. Since this is a linear least squares fit, 

we expect to find a direct solution. All that remains are some linear algebra details unique to the basis set chosen. In 
the interest of brevity, we will focus to the issue of zonal reconstruction. 

4. ZONAL RECONSTRUCTION 

We now consider the problem of zonal reconstruction for the Shack-Hartmann sensor. Here too we start with equation 
3 and derive the merit fiinction A (which is without the a,- terras): 



N 



(22) 



There is a problem with the antiderivative terms. They are quite unwelcome for a zonal fit because it would require 
discrete differentiation to recover to wavefront values. It is better to recast A to explicitly compute the wavefiront val- 
ues. So we turn our attention to the task of eliminating the antiderivatives. We begin with a close examination of 
equation 1. Note that in one dimension, the average value of the derivative of a fimction g(x) over the linear domain 
6 [a, b\ is trivial to compute: 



a ^ X _ g(^^2)-g(^>^2) 



(23) 



The trouble arises when we average in the orthogonal direction. This is what creates the troublesome antiderivatives. 
So in lieu of an exact average over an area, we use an approximation by computing the average of equation 23 evalu- 
ated at the endpoints of the domain e [c, d\ . This process is 

^x^i^vH) - 2(S3^[^^^^' ^2) "^(«» ^2))]^^ ^ ^ + (^(*. ^2) -Sici. ^2)\ ^ J 



Defining the vertex interaction matrix as 



ri = 



-100 


0 


0 1 0 


0 


0 0 1 


0 


0 0 0 


-1 



(24) 



(25) 



and the lenslet width vector h as 
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h = 



b-a 
d-c 



allows us to write the approximation as 



where the vertex vector Vis given by 



= (ia,c),ib,c),{b,d).{a,d)). 



Similar machinations yield the result 



where the interaction matrix is 



(26) 

(27) 

(28) 
(29) 



r2 = 



-10 0 0 
0-100 
0 0 10 
0 0 0 1 



(30) 



The approximations in equations 27 and 29 are now the surrogate values for the exact form of the averaged slopes. 
The appendix explores the validity of these approximations and shows that it many cases it is exact. 
At last we can replace equation 22 and pose the new merit function as 



2 

^ = W = 1 



(31) 



liberating us from the problems of having to deal with antiderivatives. The minimization is routine. The quantities 
bemg varied are the wavefront values at each vertex. All the derivatives with respect to the wavefront values at the 
vertices are set equal to zero: 



3 A = 0. 

The vector p now indexes the vertices, not the lenslets. 
The net effect is a series of linear equations that look like this (after the simplification />i » A2 = h): 

4V„,v-V„_i.v-i-V„+i.v-i-V„H.i.v+i-V„_i.v+i 

h, T 



(32) 



J^^l ^(u, v) "*" *2^(u + 1. v) ~ •*2^(a, v+l)~*l^(«+l.v+l)) 



(33) 



where f is an alignment spinor given by 



5 = 









ll 




ij 




ll 




-ij 



(34) 
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The vertex interaction matrix for this is 



^1 0 0 0 

0 ^2 ^ 0 

0 0 ^2 ^ 
0 0 0 ^, 



(35) 



where 0 has the correct dimension, i.e. 



0 = 



(36) 



Equation 33 is quite revealing in the way it relates the wavefront values Y to the focal spot shifts 8. First, the vertex 
action is between next-to-nearest neighbors, not nearest neighbors. The nearest neighbor interactions cancelled out. 
Next, the focal spot shift content is quite rich, using the measurements of the four lenslets boimdby the next-to-near- 
est neighbors. Finally, these equations are linear, so we expect a direct solution without iteration. 

This formula describes the special case for a lenslet on the interior of the array. These lenslets all have next-to-near- 
est neighbors. But many lenslets are on the edge and some are at the comers. The question is how should special 
shapes be handled? The answer is simple and comes from considering the symmetry of the problem. When a photon 
hits a lenslet, the photon can't determine the pedigree of the lenslet. The photon is deflected and impacts the CCD 
array. The amount of the deflection is the same regardless of where the lenslet is located within the array. So we are 
guided to look for a solution which respects this symmetry. 

The answer is to include an existence fimction ^ which tells us whether or not a lenslet exists. This also allows 
the user to construct arbitrary apertures. If the lenslet (u, v) exists, then ^ = 1 , otherwise e^^ ^ = 0 . Now equa- 
tion 33 can be generaUzed to solve for arbitrary masks and apertures as sculpted by the existence Wction: 

(^«, V + 1, V + 1, V+ 1 V + O^U, V 

= /'^lV.v)%v)'^'^2V+l.v)^(«+l.v)-4V.v+l)^(«.v^ (37) 



This is a series of + 1 ) (t| + 1 ) linear equations which can be written in matrix form as 

i4 • V = D. 



(38) 



Now A is the vertex interaction matrix which describes the active vertices, \\f is the vector containing the wavefront 
values at every vertex, andZ) is a vector containing combinations of the measured focal spot shifts. 

The matrix A is singular. As discussed in [4], p. 53, A is singular if there is some nullspace x that maps A to zero. 
That is there exists some vector x such that i4-x = 0. A look at equation 37 shows that each of the rows always sums to 
zero. Hence all constant vectors k are a nullspace for A. This means that we cannot find some inverse matrix A'^ such 
that 



A^'A== I. 



(39) 



But this does not mean that we cannot solve equation 38 for We are simply forced to compute a pseudoinverse ^4"*". 
For linear least squares fit problems, the method of choice to form the pseudoinverse is the singular value decomposi- 
tion, SVD. 

The process involves decomposing A into three matrices: a column orthogonal matrix U, an orthogonal matrix V 
and a diagonal matrix E. The diagonal elements of E are the singular values of i4. The process of SVD solves for these 
three matrices such that 



A = C/'I- F^. 



(40) 
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The pseudoinverse is now given by 

a"" = F.r^^/. (41) 

The entire process of reconstruction distills down to the matrix multiplication, 

yV^A'Q. (42) 

This is a very powerful statement for it shows how most of the computation can be done in advance. For a fixed lens- 
let geometry, the matrix can be precomputed and stored. It does not need to be created during the measurement 
process. 

There is some fine print however. First of all, A'^ is very large. For a 64 x 64 lenslet array, A'^ has 65"* elements 
which is about 136 MB of data. Next, the computation may take hours. So it behooves the user find ways to avoid 
recomputation when the aperture changes. Some strategies involve storing A'^ for many common apertures such as a 
series of circular masks. Another strategy is to pad the data. Consider .4'*" for a full aperture. All lenslets are active. 
Then for some reason a cluster of three lenslets is excluded jfrom the data. A good idea is to use the same A^ matrix 
and just pad the data in the excluded regions with some form of average data from the neighbors and exclude these 
vertices in the final result. 

So in closing, zonal reconstruction starts with the antiderivatives and uses matrix multiplication to compute the 
wavefront values at each vertex in one step. 



5. SOUTHWELL REDUX 

In a very influential paper published in 1980 [5], W. Southwell set forth what has become the standard method for 
wavefront reconstruction. The paper very clearly describes his reconstructor and is worth reading. Surprisingly 
though, Southwell solved his linear equations using successive overrelaxation (SOR). While SOR is an extremely 
valuable method, we feel that SVD is usually the better choice. 

We will compare the philosophical differences between the two reconstructors in the following section. Here we 
will compare the mathematics. If we define a neighbor spinor a as 



a = 



(43) 



which is just a collection of the basis vectors, then the Southwell merit function takes the form 

2 N 



(44) 



where we are now dealing only with lenslets, not vertices. In this implementation, Southwell places the computed 
wavefront value in the center of the lenslet, not in the lower left- or upper right-hand comer as expected. 

The explanation for the placement of the wavefront values is based on Southweirs treatment of the focal spot shift 
as relating to the wavefront slope - not the average slope. 

In the continuum, derivatives are measured at a point and are defined through the limit 



d^g{x^,X2) = lim 



(45) 



In the discrete space of the lenslet hy never reaches zero. It it limited by the size of the lenslet. In this space, the deriv- 
ative is defined using the difference of neighboring^points. Extending the discrete derivative is trivial. For the lenslet 
with a size vector given by equation 26, the partial derivative with respect to both variables is written 
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^1^2 ' 



(46) 



In light of this equation 3 can be viewed as 




ments, then the wavelfront values should be placed in the lower left-hand comer of the lenslet or the upper right-hand 
comer depending on whether the limit is approaching zero from the right or the left. 

Returning to the Southwell minimization problem, we proceed as before. The parameters being computed are the 
wavefront values at the center of the lenslet and we vary A with respect to these values and set the derivative equal to 
zero: 



(The existence function has been eliminated to facilitate comparison). So the Southwell interaction involves nearest- 
neighbor lenslets and uses only four measured values. 

As before we can write a matrix expression just like equation 42. The lenslet interaction matrix^ here is singular 
and we again resort to SVD. This allows for a direct solution and in many cases, we may precompute A*", Both meth- 
ods will execute in equivalent times. 



It is time to step away from the mathematical details and talk about the differences between the two reconstruction 
strategies. The new method strictly adheres to the assumptions of Shack-Hartmann operation. It treats the measure- 
ments as being extended quantities, not point measurements. This is why the wavefront value cannot be placed at a 
single point. In this model, we are only permitted to make statements about the values at the vertices. Subsequent 
mathematics are then dictated. 

In Southweirs venerable model, he assumes that the difference between the wavefront values in adjacent lenslets is 
related to the average slope across those two lenslets. One thing that seemed unusual was that the average of the two 
measurements was not placed on the boundary between the two lenslets. Regardless, the Southwell model looks at a 
lenslet, it averages the by measurement with the lenslet in the 1 direction and averages the 82 measurement with the 
lenslet in the 2 direction. 

Southwell's model mixes information from lenslets immediately. The new model assumes that the measurement of 
a focal spot shift reveals information about the differences in wavefront values at the vertices. However, the neighbor- 
ing lenslets each share two vertices. So the competing values for the shared vertices are adjudicated by the least 
squares fit. So clearly, the new model is also sharing information between lenslets, but not until the least squares fit. 

Perhaps the most obvious shortcoming with an immediate averaging is that small features of the order of a lenslet 
will be diminished. In day to day practice, this will probably not be noticed by most users. Only those requiring 
exquisite detail resolution would be driven towards the new method. 

Another interesting difference between the two is in the wavefront value interaction. The Southwell equations bal- 
ance the wavefront values with the four nearest neighbors. In the new method, the nearest neighbors interactions can- 
cel out, and the surviving interaction is between the next-to-nearest neighbor vertices. This is shown schematically in 
figure 4 below. So it is not that the nearest neighbor interaction is avoided, it is just that they exactly cancel out. 

Finally we note that the measurement content of the © matrix is richer in the new model. Southwell blends together 
four measurements, and the new model blends eight. 

In summary, both methods have a comparable run time and the memory requirements are virtually identical. They 
are both direct methods and both allow for arbitrary aperture shapes. The primary difference seems to be an improve- 
ment in resolution of fine details. Comparisons in the way the two models handle noise need to be done yet 




(48) 



For a ^ X Ti lenslet array with h^ = h2 = h, this leads to N i,r\ linear equations of the form 



6. PHILOSOPHICAL DIFFERENCES 
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Figure 4: Comp^iaglhe interactions of the wavefront values in the least squares fit. In theSotitiu^ell model depicted in (a), 
the wavefrpi^<alue is balanced between a central lenslet (marked in red) and the four nearest neighbor§>?liiscan be seen in 
equ3tki!i49. In the new model shown in (b)^ the wavefront value is balanced between a central vertex and thefeliMigxt -to- 
finest neighbors and the nearest neighbor interactions exactly cancel. The resultant interaction is specified in equation : 

7. COMPARISON OF THE TWO METHODS 
This section evaluates some sample wavefront that we put into both models and give the reader a visual comparison 
of the two strategies. The input wavefronts are shown in figure 5 below. The fust is a collection of three Gaussian 
bumps. In practice, reconstructors have trouble picking out small features against a background of big features. This 
set of bumps was designed to probe this problem. The second wavefront resembles an egg crate and it was created to 
defy reconstruction and examine the behavior of the reconstructor as it fails. 

Turning to the first case in figure 5 (a), the goal was to challenge the reconstructor in an attempt to resolve smaller 
bumps in the presence of a large bump. Here we see two rather sharp bumps in firont of a prominent feature. The 
smaller bump on the left is taller and wider than the bump to the inmiediate right. This bump is the most difficult one 
to resolve. 

The results are shown in figure 6 (a) in the top row. Since the smaller bumps are the features of interest, the plots 
have been scaled to their height. The input wavefront has been plotted at the same resolution as the lenslet array to 
facilitate connparison witih the reconstructor outputs. As expected, the new reconstructor fared better at resolving the 
small bumps. The bump on the left is faithfully reproduced. The bump on the right is reproduced vaguely. All recon- 
structors have an amplitude ambiguity, so we quantified the problem by essentially doing a least squares fit to see 
what magnification it would take to get the reconstruction to most closely match the input. The new reconstructor 
overstated the amplitude by 0,4%; the Southwell reconstractor understates the amphtude by 2.7%. The values 
were 0.003 for the new method and 1 1 .6 for the Southwell method. 




Figureor: iwo omerent wavetronts are shown here to highlight the differences between the two reconstructors. TheT»fet case 
t^in (a) presents three bumps with different sizes. This shows how the. two methods are able to resolve features of different 
sizes. The second case in (b) was created to be difficult or impossible to reconstmct. The results were surprising. 
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N^O 30 40 

new reconstructor Southwell rfesQnstructor 

Figure 6: Results comparing the new reconstructor to the Southwell reconstructor. The wavefronts shown in figure 5S^iere 
submitted to the two reconstructors. The top row (a) shows the input wavefront at the resolution of the lenslet array and 
shows the results of the two reconstructors. The new method recreates the smaller bumps better indicating that this method is 
best suited to resolving fine details. The bottom row (b) shows a a series of two dimensional plots because at the resolution of 
the lenslets, the three dimensional plots are difficult to interpret. Surprisingly, the output from the new reconstructor matches 
very well with the input. The Southwell method did note fare as well, but still performs adequately. This suggests that the new 
reconstructor is better suited to handled complicated structure. 

After running dozens of cases, we wanted to present a wavefront that could not be reconstructed in order to study 
the behavior as the algorithm failed. We wanted to identify artifacts associated with an errant reconstruction. We cre- 
ated the wavefront shown in figure 5 (b). The magnification computation shows that the new reconstructor overstates 
the amplitude by 1 1% and that the Southwell method imderstates the amplitude by 19%, The values were shock- 
ing. The new method has a of 2 10"^^; this represents the limit of double precision rational. In other words, the new 
method had a %^ of zero. For the Southwell method the value is 11,800. 

The fact that the new reconstructor had essentially no fitting errors was astounding and led to a careful exploration 
of the approximations in equation 27 and 29. The results, detailed in the appendix, show that the approximation 
allows for a very rich wavefront structure across each lenslet. It probably allows for more structure that the focal spot 
location method can tolerate. That is, the focal spot will degrade and become ambiguous before the approximation 
used in the reconstructor breaks down. 



8. SUMMARY 

This paper details new reconstruction strategies formally based upon the principles of the Shack-Hartmann sensor. 
This- method was compared to a modernized implementation of the standard method developed by Southwell. The 
difference between the two different procedures is negligible or extreme depending upon the application. For com- 
mon measurements of smoothly varying wavefronts, the difference between the two methods teeters on the edge of 
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imperceptible. However if the wavefront of interest is richly structured or finely detailed, then the new method is a 
substantial improvement. 

APPENDIX: VALIDATING THE APPROXIMATION 

The results on the egg crate reconstruction (figure 5 (b)) were startling. After all, the wavefront was designed to break 
the reconstnictor, yet the was zero to the limit of machine precision. This inspired a detailed look at the approxi- 
mation used in equations 27 and 29. This section documents the investigation. 

To recap, the merit function (equation 22) for Shack-Hartmann sensors involves antiderivatives, which are unwel- 
come in a zonal reconstruction. An approximation was used to bypass the nettlesome antiderivatives and formulate a 
merit function based on the wavefront values. The approximation is contained in the two equations 

i-J^^^^-H^W^- (50) 

We start by forming the sum and the difference of the two equations represented above. Remembering that the 
domain is Xj e [a, i] , e [c, d] we find that 

\|/(6,^f)-V|/(a.c) = f □^2(^1,^2) + fn4'i(;ci,jC2). (51) 

and 

\(f(i^,c)-M/(a, J) = jj-a^2(^i»^2)~r-°'^'i(^i'^2)- (52) 
"2 "1 

At this point we need to be more specific about ^(xp ^2) . We can represent this function as a Taylor series. This is 
not because we champion the use of the Taylor monomials in modal reconstruction; it is because any piecewise con- 
tinuous function can be described exactly by an infmite series of the monomials. We now write 

When this expansion is substituted into equations 51 and 52, we see 

d / 

i-oy = o 



i-0; = 0 i = Qj~0 



(54) 



and 

d i 
i « oy « 0 



d i 

a. 



i = oy = o i-oy»o 
Using the identity 
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^:rt X" ^ (56) 

one can show that the approximations in equation 50 are exact whenever the wavefront over the lenslet has a structure 
like 

oo oo 

^2) = «00 + «ll^l^2 + X «'0^'l + E «0.-^2 • (57) 

This allows for an incredfcle amount of stnicture over the lenslet. For example, any separable function 

VC^i.Xj) =^(x,)F(x2) (58) 

will be hancfled exactly; in other words, there is an exact solution and it is described by equation 42. Or any separable 
Junction with the lowest order torsion term, for example 

Xj) = ^(ac,)7{x2) + kx^x^ . (59) 

The egg crate wavefront in figure S(b) is 

VC^l.Xj) = sin(A:iXi)cos(A:2A:2) (60) 
which is exactiy the form of equation 58. 
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[000186] EFFICIENT COMPUTATION WIH SPECIAL FUNCTIONS LIKE THE CIRCLE 
POLYNOMIALS OF ZERNIKE. 

[000187] The circle polynomials of Zemike play a prominent role in optical analysis. While 
decompositions of wavefironts into Zemike polynomial series can yield valuable insight, 
computing with the polynomials themselves is quite inefficient. Here we outline how rational 
polynomials like those of Zemike, Legendre, Chebyshev and Laguerre can be handled as afSne 
combinations of a Taylor monomial set. We demonstrate how calculations can be performed 
much more rapidly in the Taylor basis and how to use integer transformations to recover the 
exact amplitudes in the desired basis. We also explore CH-+ optimizations for storing the Zemike 
amplitudes and transforming between Zemike polynomials and Taylor monomials. 
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1. INTROIJDCnON 

Urn dide pdjuamials of Nobel laureate Fritz Zenuke have proven to be an invaliMble tool far investigatdiLg optical 
aystBms* Zeniks's pdynoioialfi are csrthoganal which, fkiilitaties assodating individual tesnm vritii physical obaerv- 
ables. GenaratiQiis of ojitical ssieaitists haveused tihis represenialiaEito desGrihQ optical system pafbmiaiice, ta quan- 
tify behavior of individual compoaenls and to provide a fiamfiwcxrk for conasiEat ovaliaatioa of oprtdcal phmomianAa 
All ftiB qtianfeificafioa is done in terms of flie amplitudes fe eadi teim in a Zemike polynomial expansioiL While 
prnike'te polyaaruials have proven to be an invaluablQ medianisni for describing the physical warld^ fee feomdatiQn 
is qnite cxunbersame and actual camputatidas in this basds set are BJewdlessly camplicated. 

Becausa^of the Qram-Sdanidt arthpeonalizationprocess feat ZemikB used to construct his pdynornials, fee ccsm- 
pleodty of higher ardear terras grows rapidly Fortunately^ the Taylor monomials form a minimal complexity rap- 
resentafioa that is exactly equivaknt ta tie ZiemikB. basis on an oidfir-by-arder basis. This paper fbnnally estabUshfis 
ftifl equivalence of fee two sets and ten discusses how to con^sutB in tiie l^ykar basis and transform the results into 
the uaafiil and femihar ZHrnika basis. In fact^ this procedure is evenmore general: the Taylor mmoniials can be used 
to r?presQEA any rdaaual pdynotnial se** sudi as feosa of Legendrej Laguen^e and CJheljj'shfiv. With ftiis realization, 
ftie Taylor basis becomes ev^ raare valuable for it is hi essence all flis otier rational pdynmdal basis seta. So by 
coniputdnginthfi ISiylor hasia, ane is ablato useaiSne.transformations to recoverlheaiiiplitudes of aU other rational 
pdyncinial basis. This nieans, fbr example, "thatthere needs to he only one recqnstructor Wlfti the Taylor amplitudes 
in. hani we are^ only caie matrix mtaltfplicatian away from the amplitudes of Zemike^ liegendie, Oiebysliey oi' 
Lagume. And given the arnpli^tucjes in tiaasB other satej a matiism ten into Taylor amplitudfis, 

and feea anofiier maMxmdtipiiGatipn can tdca us iiato the basLyset of dboice. 

TMsjpaper discuss some praddt^l^exa^ Zermk© polynomials aracurnbei'some 

and tedious to use. We explain how to use the Ujlor basis for cornput^an and then how to recovo: fee Zeniike 
amplitudea. Besides laying out fe& Jtodarniantal maftieniaticiB., we go into to detail to show ho w Ms would be imple- 
meated in, a modm cdmpUtBr lEftiguiage sudh as OH-, Forefhouight in diBsi|ping fee aiciiifasduro of code results in 
an extremely fadle andpowerful tool fbr wayefrontaiialy.sLs. 

Thematheimfesr of abstractamjandliia CQmptiter code shows a djefiiutB^x^^ 

eral issues in implementatdmi sujeh as run time aud memoa^y usage,, are explored. Finally tasHng and debugging a 
laigp volume ofcodei-equires stringent validation. We outhne tiie five layers of validation we used to verify several 
znil liaa lines of computer code;. 

[000188] 
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2. EQUIVALENCE OF THE ZERNIKE AND TAYLOR AMPLITUDES 

When the complete set of Zemike^ and Taylor functions through arbitrary order d are considered, the two sets are 
equivalent. Consider the explicit example of two sets for second order. 



Taylor function 


Zernike function 


^oof^^ y) = 1 


^oo(^' y) = 1 




Zio(x. yj^'x 


Toi(x. y)=y 


'Z'\\(x,y)^y 


hoCx.yJ'^x^ 


Z2o(x.y) = 2jcy 


'^n(x,y)^xy 




To2(x>y)-/ 


Z22(^. y)-y'-x' 



The equivalence is shown below. First, we express the Zernike functions in terms of the Taylor functions: 

Zoo(x. y) = Tqq(x, y) 
2io(x, y) ^ Tio(x, y) 

Zn(x,y) = T^x(x,y) 
'Z"i^(x,y)-rr^x(x.y) 

221 (x. y) = 2T2o(x, y) + 2T^2(x, y) - Too(x. y) 

y) = 702^^^ y) - T2(i(x. y) 

Next we express the Taylor functions in terms of the Zernike functions: 

Toq(x, y) = Zoofx, y) 
'^\o(X' y) = Ziq(x, y) 
T^^\(x,y)^1\\(x,y) 

T2^(x. y) = -\ Z22(x, y) + \ Z2ifr y) + \ Zoo(x, y) 
Tn(x.y)^\z2^(x,y) 

To2(x. y) = I Z22r^, y)'^\z2i(x.y)'^\ z^(x^ y) 

So for example the arbitrary wavefront , written as 

y) = too + t^QX-^ t^^y + t^QX^ + t^^xy-^ t^^y^ (1) 
can be expressed in terms of the Zernike polynomials as 

- ^^.0-.3'<^^4>,x.,-H(^^-^f)z,.., (2) 

One already sees an inkling of the native simplicity of the Taylor basis. The relationship between the two bases can be 
cast in matrix form as 



b. The original notation of Zernike has been reproduced in the classic work by M. Bom and E. Wolf, Principles of 
Optics 7e, Cambridge, 1999. However, we find the notation of Malacara [2] to be more practical here. 
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} 



1 0 0 0 0 0 
0 1 0 0 0 0 
0 0 1 0 0 0 
0 0 0 0 2 0 
-1 0 0 2 0 2 
0 0 0 -1 0 1 



^00 

^20(^'3') 

y) 



(3) 



^00 (^»>') 
T^o^^'^y) 
Toi{x,y) 

T2oix^y) 

Tn{x,y) 
To2ix^y) 



10 0 0 0 0 

0 10 0 0 0 

0 0 10 0 0 

1 0 0 0 i 

4 4 2 

0 0 0^0 0 

1 0 0 0 i \ 

4 4 2 



^Qo^^^y^ 
^ioi^.y) 

Zii{x,y) 

ho^^^y) 
Z2[{x>y) 

Z22(X,7) 



(4) 



But even this notation is cumbersome and it is more convenient to write these relationships in the forms 

Z(x,y) ^fz^.T{x,y), 



(5) 



T{x,y) = tz'^Z{x,y). 



(6) 



The matrix tz transforms an amplitude vector from the Taylor basis to the Zemike basis and the matrix /z transforms 
an amplitude vector from the Zemike basis into the Taylor basis. Their usage in equations 3 and 4 may seem counter- 
intuitive, so a short demonstration follows. 

The mappings in equations 3 and 4 show that the Zemike and Taylor bases are equivalent. This means that any 
wavefront y\f{x,y) can be resolved in either basis allowing us to write either 



\\^{x,y) = t T{x,y) 



(7) 



or 



y]f{x,y) = z •Zix.y) 



(8) 



where t and z are respectively, amplitudes in the Taylor and Zemike spaces. These amplitudes are the outputs from 
modal reconstructipn. Notice that the results of the dot products are a scalar function, in fact, ttiey are exactly the 
same scalar function as seen in equation 2. Since equations 7 and 8 are the same function, we can write 



z^'Z(x,y) = t^'T{x,y) 



which in light of equation 5 can be recast as 



This immediately implies 



which is of course 



z'^'ifz^nx.y)) = t'^nx.y). 



T -T T 
Z fz = t 



t =fzz. 



(9) 
(10) 

(11) 
(12) 
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A similar process yields the relationship 

z =^ tz't. (13) . 

The naming scheme for the transformation matrices should now be clear. The matrix tz takes amplitudes from the 
Taylor basis to the Zemike basis and fz takes amplitudes from the Zemike basis into the Taylor basis. As one might 
expect, these operations are the inverse of one and other. That is 

tz 'fz ^fz-tz=l (14) 

where 1 is the identity matrix of appropriate order. 

This analysis can be extended arbitrarily order by order. The approximation theory of Weierstrass [5], p. 1929 
shows that in the limit d-^oo either set of polynomials can exactly represent any piece wise continuous wavefront. 

The matrix^z is based upon the recursion relationship in [2], p. 471: 

q m m—j 

1 (x,y)= Y Y Y (~1)''^Y""^^¥^^^-^1 ^2(/ + *)+p n-2a+y+A)-p . . 



where s = sgn(v), p = '|^(j+l),i7 = |(v-~ 5Mod2n) , and = ^ - m - |j . For a given order n and index m, this 

formula can be used to generate a vector of coefficients for the fz matrix. When the fz matrix is ftilly assembled, the tz 
matrix is computed as 



tz = fz 



-1 



(16) 



Of course the transformations represented,^ and te, are not just a single matrix; they are a series of matrices, one pair 
for each order. 

By comparison, the Taylor monomials are generated by a much simpler function: 

T^^jj{x,y) - x''-^}/ , (17) 

In summary, it is perhaps convenient to consider the two bases, Taylor and Zemike, as having a relationship simi- 
lar to the one between the Cartesian and spherical polar coordinates. Certainly we all realize that (x, y, z) references 
the exact same place as (r, 0, <|)j. These mappings are both complete and unique. They are completely equivalent ways 
to express location. Integrations in both coordinates yield identical results. 



3. THE FUNDAMENTAL NATURE OF THE TAYLOR MONOML^LS 
We are all familiar with the Taylor expansion and we know that any piecewise continuous ftinction can be exactly 
reproduced by an infinite number of expansion terms. The Taylor monomials represent the simplest basis for an 
expansion. All other expansion sets, such as those of Laguerre or Legendre, are just combinations of the Taylor 
monomials. The natural question then is why are these other bases so popular? It is primarily due to the fact that the 
expansion polynomials are orthogonal. In general, two functions ffx^y) and fj(x,y) are orthogonal over some domain 
D with respect to a weighting function w(x,y) if they satisfy 

l\fiixyy)fj{x,y)w{x,y)dxdy = a^jS.j (18) 

D 

where ay is some constant and 6^^ is the Kronecker delta tensor. 

We are all familiar with the concept of orthogonality since we live in a Euclidean world. Consider some point P on 
an x-y graph. We are able to measure x andy independently because the axes are orthogonal. As we vary in the x- 
direction at P, the y values do not change. In general for a coordinate system we can formulate the geodesic between 
two points by integrating ds^ using the metric g 
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) 



'J 



as shown in ref [1], p. 105. For the Euclidean world that we live in, dqi = dx, dq^ = dy, and the 



metric is 



1 0 
0 1 



(19) 



(20) 



hUhis special case, and rfy may have finite size. In general, orthogonal coordinate systems have gy = 0 when i 

SiJ = ^iAj (21) 

where the by are constants or functions.This states mathematically what our intuition tells us: we can be at a point P 
and vary our position m one coordinate without affecting the values in the other coordinates 

A companson of equations 1 8 and 20 suggests the concept of orthogonality extends in a natural way to functions 
The consequence IS tiiat measurements of the amplitudes are independent of each other. For example, consider the 
Zemike polynomial. It is possible to do a least squares fit to measure a single amplitiide like Z42. The measurements 
are mdependent and a computation of subsequent amplitudes will not affect tiie previously measured amplitudes On 
flie otiier hand, a non-orthogonal set of basis function like the Taylor monomials behaves differently The amplitiide 
for a term like^ / depends upon which tenns are being measured with it. This makes the matter of physical interpie- 
^tion quite difficult. Consider the lowest order Zemike polynomial term. Zoo- It's amplitude is the mean value of the 
^^eS^frSS fit "^^y^"' amplitude - the constant term - will fluctuate depending upon how 

nrth^r."?! ' t' '! '° f ° computations in a tiie Taylor bases and then transform the amplitudes into an 
orthogonal set. The Taylor monomials are the simplest bases set to work with. Figure 1 shows schematically how the 
fimctioi^l bases are related. The Taylor basis is the hub and is the set of minimal complexity. After doing computa- 
tions m the Taylor basis, rational transformations are used to compute quantities in any desired basis. This also allows 

^nZl^T'''''' ''T'L'''^ T "^'^ ^ Chebyshev (with the usual caveat that the domains 

must match) For example, if/c and cz are the transformations that take amplitiides from and to the Chebyshev basis 
then we could start with a vector of Zemike amplitudes, z, and compute the vector of Chebyshev amplitiides, c, via ' 

(22) 



c = tcfzz. 



For the case of « polynomial bases, there are 2n tiransformations. If we don't use the Taylor bases as tiie hub. and have 
tiransfomiahons direcrty comiecting all the bases, then n(n - 1 )/2 transformations are needed. The Taylor hub will 
greatly sunplify tiie affme tiansformations. 




Figure 1: Think of the Taylor monomials as a hub comiecting all the rational polynomial basis sets; four such sets are shown 
here.W.th software deigned around the Taylor basis, the user assures the fastest possible execution and the least amount of 
computer code. An affine toolkit is used to transform the resulte of the Tkylor-based calculations into tiie other basis sets 
These atfine transformations are integer based and virtually instantaneous. 
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(a) (b) 

Figure 2: Aperture resizing is just rescaling the coordinate system. Consider an arbitrary function shown in (a). We can shrinlc 
the ;c-coordinate by a factor of/= 1/2, as shown in (b). The curve does not change. We have efifectively just relabelled the 
coordinate system. 



i = Oy = 0 



(26) 



Of course practical applications truncate tiie e^qjansion at some finite d. Now we apply the transformation in equation 
25 to equation 26 and we see that the amplitudes of order d are multiplied by the factor / . For example the transfor- 
mation matrix F that adjusts the amplitude vector for a second order expansion is 



F = 



/ 0 0 0 0 0 
o/ 0 0 0 0 
0 O/ 0 0 0 
0 0 O/ 0 0 
0 0 0 o/ 0 
0 0 0 0 o/ 



(27) 



Now if the user had Zemike amplitudes foisted upon him and wished to perform the resize operation, the affine tool- 
kit will help. Given a Zemike amplitude vector 2 and a resize matrix F, the new Zemike amplitudes z ' would be given 
by 

z' = tzFfz'Z. (28) 

And were we given Zemike amplitudes and asked to reduce the aperture and calculate the Chebyshev amplitude vec- 
tor c we would compute 

c = tC'Ffz'Z, (29) 

Yes, it may be possible to compute F in other bases, but it is trivial in the Taylor basis. And working in this basis 
allows the user to compute amplitudes in any other rational polynomial basis set. 

5. SPECIAL CASES 

While in general the Taylor polynomials are the most efficient computation basis, the rule is not absolute. There are 
instances when the evaluation is much faster in the Zemike basis. These cases involve wavefronts with rotational 
symmetry, such as spherical waves. 
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To study how increasing the fit order d improved the resolution of sharp details, a cone was modelled. The equa- 
tion for a cone in polar coordinates is 



Y(r,G) = 1-r. 



(30) 



The strategy is to perform a least squares fit to find the amplitudes which best describe the cone. If we compute the 
Zemike amplitudes, then we may compute them one order at a time. The orthogonality of the Zemike polynomials 
guarantees that as higher order amplitudes are computed, the lower order amplitudes will not change. This is the great 
advantage of an orthogonal basis set. The merit fimction for order n is written as 



0 0 



m = 0 



rdrdQ . 



(31) 



The merit fimction is minimized by setting the derivatives with respect to the amplitudes equal to zero and solving the 
resultant set of linear equations. For example, picking the arbitrary amplitude c,y, 



= 0. 



This immediately leads to the equation 

2jc1/ n 



0 OVm-O 



27Cl 



Zy(r, Q)rdrdQ = J j\i/(r, e)Zy(r, Q)rdrdd . 



(32) 



(33) 



0 0 



There is very little integration needed to solve this problem. The left hand side is solved by using the orthogonality 
relationship (a corrected form of the equation 13.11 in [2]) 



(34) 



and we see that this integral projects out the value of c,y. The right hand side is greatly simplified because of the peri- 
odic nature of trigonometric fiinctions. Since sin(9 + 2n) = sin(0) and similarly for cos(e), the only integrals that 
can be non-zero are the rotationally invariant terms which do not have a 6 dependence. The general solution for the 
only surviving coefficient is given by 



c n 



2nl 

n + l 



^Jjv(r,e)Z ^{r,Q)rdrdQ. (35) 



^ 0 0 2 

This greatly simplifies the problem. In a few minutes we were able to perform a a Zemike polynomial fit to order 
1000. This expansion series looks like 

2 2 2 2 

437^20, 10 " 525^20. 10"*" 10197^100. 50" ••• 1001997^1000.500 • (36) 

Attempting this feat in the Taylor basis would require the inversion of a 501,501 x 501,501 matrix, so this method 
was much faster. The issue here is that this problem is analytic, not numeric. Applications involving measurements 
are typically numerical, and in those cases the Taylor basis is much faster and easier. 

Figure 3 shows some results for this fit. The first image shows the input fimction in equation 30. Subsequent 
images show the fits to successively higher orders. The first order term, the piston, is the mean value of the function. 
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The higher order terms add curvature. Eveatually they will sharpen the point and straighten the sides. By order 40, the 
fit closely resembles the input curve. 

At this point, we hope to have estabhshed the background for and strongly motivated the use of the Taylor mono- 
mial basis set. The final component is to discuss the issues of putting these ideas into an efficient, rapid and robust 
computer code. Some forethought is required, but after starting with a good design, the code is very easy to use. 

6. IMPLEMENTATION 

6.1 Software Design Issues 

To support the aforementioned methodologies of computation, data structures and algorithms need to be optimized 
for this solution. Because of the flexibility and speed requirements of these computations this implementation is done 
in C++. 

A three dimensional array can be used for the representation of the Zemike polynomials. By using a three dimen- 
sional array, it becomes simple to access any defined order, polynomial or coefficient. The first dimension of the 
array is the order and is addressed Z[o], The second dimension of the array is the polynomial and is addressed 
Z[o][p]. Finally the third dimension of the array is a coefficient, which is addressed Z[o][p][c]. 

Usmg a three dimensional array can be very inefficient if not implemented correctly. As orders increase, so do the 
number of polynomials and coefficients. To maintain our addressing structure, a statically allocated array would be 
very wasteful. Consider building the Zemike array to third order. A statically allocated array would have to be 
declared INT64 Z[3][3][10]. This array statically allocates 3 polynomials and 10 coefficients for all orders even 
though orders 0, 1, and 2 have fewer polynomials and coefficients. By using dynamically allocated memory we can 
prevent ourselves firom using excess memory and still maintain our addressing structure. Building the Zemike stmc- 
ture up to but not including nOrder is actually very trivial. Listing I below shows a sample allocation. 




Figure 3: An example of using the Zemike basis to speed computation. Starting with a rotationally invariant object, the cone 
shown in the upper left-hand comer, successive least squares fits were done to extract the amplitudes of the rotationally 
invariant Zemike polynomials. A fit to order 1000 took only a few minutes. This computation may not be feasible using the 
Taylor basis. The results of fits to various order are shown in subsequent frames. 
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'typedef INT64 Coeff; 
typedef Coeff* Poly; 
typedef Poly* Order; 
typedef Order* zernikes; 

int nPoly; 
int nCoeffs; 

try{ 

2 = new Order [nOrder] ; 
for( int i=0; i<nOrder; 
nPoly = i + 1; 
nCoeffs = ( (i+l) * (i+2) 
z[i] = new Poly[nPoly] ; 
for (int j=0; j<nPoly; j++ 

z[i] [j] = new Coef f [nCoeffs] ; 
memset(z[i] [j] , 0, sizeof (Coef f) * (i+l) ) ; 

} 

} catch ( std: :bad_alloc ) { 
FreeZernikes (z) ; 
return ; 



i++ ) 



{ 

) / 

) { 



2; 



Listing i 



// Coeffs 64-bit integers 
// Polys array of Coeffs 
// Orders array of Polys 
// zernikes array of Orders 

// Num poly 
// Num coeffs 

// Try to allocate 

// Allocate orders 

// For each order 

// Polynomial in order 

// Coeffs in polynomial 

// Allocate polynomials 



For each polynomial 
Allocate coeffs 
// Set coeffs to zero 
// Matches inner for 
Matches outer for 
Catch alloc error 
// Free memory 
// Exit function 
// Matches catch 



// 
// 



// 
// 



Constructing the Zemike coefficient table in this fashion avoids unnecessary allocation of memory. After construct- 
ing this structure, the non-zero coefficients need to be set to the correct value. This can be done in one of two ways. 
First by calculating the non-zero coefficients for every polynomial in all orders. The following code segment gener- 
ates these non-zero coefficients for any Zemike polynomial where n is the order and m is the offset of the Zemike in 
order n. 



INT64 C, p, d, s; 

double M, q; 

long mu, nu, kappa; 

d = n - 2*m; 
s = Sgn(d) ; 

M = (n/2,0) - fabs(m - (n/2.0)); 
p = (abs(s)/2.0) * (s + 1) ; 

q 



s * (d - s*(n%2) ) / 2.0; 



for (long i = 0; i <= q; i++) { 

for (long j = 0; j <= M; j++) { 

for (long k = 0; k <= M- j ; k++) { 
mu = 2 * (i+k) + p; 
nu = n - 2 * (i + j + k) - p; 
kappa = mu + nu; 

C = pow(-l, i+j) * Binomial (n-2*M, (2*i)+p) 
Binomial (M- j , k) * (nFact (n- j ) / (nFact ( j ) 
nFact(M-j) * nFact(n-M-j))) ; 
INT64 index = (kappa * (kappa+l))/2 + nu; 
Z[n] [m] [index] += c; 



} 



} 



// Index variables 

// Angular frequency 
// Side of periodic table 
// Dist from center line 
// cos/sin term determinant 
// Loop variable 



// Update mu 

// Update nu 

// Update kappa 

// Evaluate a coeff 



// Coeff position 
// Update the coeff 
// Matches inner for 
// Matches middle for 
// Matches outer for 



Listing 2 
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It*s obvious that this calculation is at least a @{n^) operation due to the multiple embedded for loops. The second 
method is to load the non-zero amplitudes into the array manually, the previous section of code can be easily modi- 
fied to generate the C-H- code to make this a 0(1 ) operation. 



z[03 [03 [0] 



z[l] [0] [1] 
z[i] [1] [2] 



(INT64) 1; 



(INT64) 1; 
(INT64) 1; 



// 0th order coeff 

// n = 0; m = 0; constcint 

// 1st order coeffs 
// n = 1, m = 0; X 
// n = 1, m « 1; y 

// 2nd order coeffs 



z[2] [0] [4] 




(INT64)2; 


II 
/ f 


n = 


2, 


m = 0; 


xy 


Ty^i r^i 

z[2] [1] [0] 




- (INT64)1; 


II 


n a 


2, 


m = 1; 


cons 


2 ^2^ rn r^i 




\ JLJM 104 } ^ t 


II 


n = 


2, 


in « 1; 


x» 


zC2] [1] [5] 




(INT64)2; 


II 


n = 


2, 


m = 1; 




z[2] [2] [3] 




' (INT64)1; 


II 


n = 


2, 


m = 2; 


x^ 


z[2] [2] [5] 




(INT64)1; 


II 


n s 


2, 


m = 2; 


y^ 








II 


3rd 


order coeffs 


z[33 [0] [6] 




-(INT64)1; 


II 


n = 


3, 


m = 0; 


x» 


z[3] [0] [8] 




(INT64)3; 


II 


n 


3, 


m = 0; 


xy* 


z[3] [1] [1] 




- (INT64)2; 


II 


n = 


3, 


m »= l^- 


x 


z[3] [1] [6] 




(INT64)3; 


II 


n = 


3, 


rn = 1; 


x» 


z[3] [1] [8] 




(INT64)3 7 


II 


n =s 


3, 


m = 1; 




z[3] [2] [2] 




- (INT64)2; 


II 


n = 


3, 


m = 2; 


y 


z[3] [2] [7] 




(INT64) 3; 


II 


n = 


3, 


m = 2 ; 


x^y 


z[3] [2] [9] 




(INT64)3; 


II 


n = 


3, 


ra = 2; 


y 


z[3] [3] [7] 




- (INT64)3; 


II 


n = 


3, 


m = 3 ; 


x*y 


z[33 [3] [9] 




(INT64)1; 


II 


n = 


3, 


m = 3 ; 





// 4th - nth order coeffs 



Listing 3 



Since this methodology uses only assignment operations, and lacks any mathematical computation or for loops, this 
segment executes in a constant time. The second method makes a larger data segment of an executable, DLL or a stat- 
ically linked library. To illustrate the differences, a sample program that builds and frees the Zemike array 1000 times 
was created and the coefficients were populated using both methods. The results are summarized in table 1. The first 
method creates an executable that is 168 KB. The same sample program using the second method creates an execut- 
able that is 2 MB. When these two programs were ran on a machine with two 1.5 MHZ processors, the first method 
took on average 385 ms while the second method only took 14 ms on average. Clearly speed is achieved by using 
more disk space. When speed is an issue the second method should be used. 

One other issue comes into play while deciding which method to implement. The coefficients that represent the 
Zemike polynomials through order 39 are small enough to be represented by a 64-bit integer. However, due to the use 
of factorials in the computation of these values, an arbitrary precision number type is necessary to compute them. In 
fact, using a naive implementation widiout arbitrary precision, the maximum order of the coefficients that can be cor- 
rectly computed is 20. To calculate the coefficients during runtime, ref [3] outlines fimctions necessary to write a C 
implementation that handles extremely large numbers and the algorithms they chose are highly efficient. Otherwise, a 
tool such as Mathematica will prove itself priceless for generation of the coefficients, since all of its calculations are 
done using arbitrarily precise numbers. It is also necessary to note that above order 40, the coefficients become large 
enough to overflow a 64-bit integer data type. 
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Method 


Time to Load 
Amplitudes 


Memory Used 
During 
Execution 


Memory Used 
on Disk 


Computing 


385 ms 


1388 KB 


168 KB 


Amplitudes 








Pre Computed 


14 ms 


3312 KB 


2MB 


Amplitudes 









Table 1: Methods for Loading Amplitudes 



Now that we have discussed the internal representation of the Zemikes, the computation transforming them into 
their Taylor representation needs to be considered. The matrix operations in equations 3 and 4 are an example of how 
the transformations are performed for order 2. Matrix multiplication quickly becomes very expensive. For example 
the transformation matrices at order 30 are 496 x 496. However, since the matrices are very sparse, it is extremely 
time and memory inefficient to create the transformation matrices and then perform a multiplication. Using a similar 
techoique as the one used to create the Zemike vectors, we use equation 15 to create code to perform sparse matrix 
multiplication. For example, converting a 4th order Zemike polynomial to its Taylor equivalent can be perfomied as 
follows: 



tay[03 




zer [0] 


+ zer[3] / 


(INT64)4 


+ zer [5] / 


(INT64)4 


tay[l] 




zer [1] 


+ zer[6] / 


(INT64)2 


+ zerES] / 


(INT64)6 


tay[23 




zer [2] 


+ zer [7] / 


(INT64)6 


+ 2er[9] / 


(INT64)2 


tay[3] 




zer [4] 


/ (INT64)2; 








tay[4] 




zer [3] 


/ (INT64)4 


+ zer [5] 


/ (INT64)4; 




tay[5] 




-zer [3] 


/ {INT64)2 


+ zer [5] 


/ (INT64)2 




tayt6] 




-zer [6] 


/ (INT64)4 


+ zer[8] 


/ (INT64)4 




tay[7] 




zer [6] 


/ (INT64)4 


+ zer [8] 


/ (INT64)12 




tay[8] 




zer [7] 


/ (INT64)12 


+ zer [9] 


/ (INT64)4 




tay[93 




-zer [73 


/ (INT64)4 


+ zer [9] 


/ (INT64)4 





Listing 4 



where tay is the array holding the Taylor coefficients and zer is the array of Zemike coefficients. With this imple- 
mentation, the transformation of the Zemike coefficients to the Taylor representations performs only 19 divisions, 
whereas the matrix form would do 225 divisions. To do the inverse all that is necessary is to perfomi these operations: 



zer [0] 




tay[0] 




tay [43 ; 




zer [1] 




tay [13 




tay [7] * 


(INT64)2; 


zer [2] 


s 


tay [23 




tay [83 * 


(INT64)2; 


zer [3] 


ES 


tay [43 


* 


(INT64)2 


- tay [5] ? 


zer [4J 




tay [3] 


* 


(INT64)2; 




zer [53 




tay [43 


* 


(INT64)2. 


+ tay [53 ; 


zer [6] 




-tay[63 




f tay [7 3 * 


(INT64)3; 


zer [7] 




tay[83 


* 


(INT64)3 


- tay [93 * (INT64)3; 


zer [83 




tay [63 


* 


(INT64)3 


+ tay [73 * (INT64)3; 


zer [9] 




tay [8] 


* 


(INT64)3 


+ tay [93 ; 



Listing S 



Once again, very few math operations are required compared to the matrix multiphcation method. Thus, it becomes 
superfluous to use large amounts of memory to hold the transformation matrices. Also, since computation of the 
matrix values isn't performed at runtime, the transformations are extremely fast, ©(1) . There is one disadvantage, 
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hard-coding the transformations can create extremely large executable files; approximately 60 MB for both transfor- 
mations through order 40. 

Plotting a wavefront using the Taylor monomials becomes trivial. First the user would define the Zemike polyno- 
mial in a vector format. Next the amplitudes would be transformed to Taylor amplitudes. Finally we can use tiie Tay- 
lor representation to plot the wavefront 

c = (0, 0, 0,0,1,0) 
a=fz*c 

n 

i" 1 

where p is the vector holding the exponents for x and y. By using a class object we can write C++ code that looks 
almost identical to the mathematics above. Since all transforms will support the same interface, many different trans- 
forms can be defined by using inheritance. If a class object is defined for transforms, we can then define what it 
means to multiply a transformation matrix by a vector of amplitudes by using operator overloading. The following 
C++ code does the same operation as the mathematics above 



CTransform *fz = new CTFZ(); // Define transform 

Wavefront W; // Define wavefront 

Poly c = { 0, 0, 0, 0, 1, 0 }; // Define Z42 

Poly a = fz*c; // Transform to Taylor 

for( int i=l; i<=n; i++ ) // sum 

W[x] [y] = a[i]*pow(x, p [i] [1] ) *pow(y, pti][2]); // Taylor components 

delete fz; // ^^^^ transform 

Listing 6 

First the code creates the class that transforms to the Taylor basis from the Zemike basis. Next we define a wavefront. 
The polynomial c is created to represent the eigenstate Z4 2 and then the Taylor amplitude vector a is created by mul- 
tiplying c by the transform. Finally the wavefront is calculated in the Taylor representation. 

6.2 Code Validation 

Using the said methods for code generation can yield hundreds of thousands of lines of code, depending on the order 
desired. Debugging that many lines of code is quite difficult, therefore, the use of mathematical validation is unavoid- 
able. 

Several tests were done to show that the Zemike amplitudes were computed correctly without a single missed 
keystroke. For these tests, an independent programmer manually wrote C++ fimctions for the first five orders of the 
Zemike polynomials. So the production code and validation code were created by two different programmers. The 
first validation test was a point by point test. This test compares the fimctional definitions of the Zemike polynomials 
through fiftii order with the results of from using the Zemike array stracture. The point by point test samples 10^"^ 
points within the unit circle and compared the values returned for both methods. There are no differences between the 
fimctional or the array representations. The next test uses the property of orthogonality to determine if the Zemike 
amplitudes were created correctly. The orthogonality of the Zemike polynomials, stated in equation 34, shows that if 
the product of two Zemike polynomials is integrated over the unit disc, the answer will be zero unless both polynomi- 
als are the same. The orthogonality test was comprised of three parts. The first part tested the orthogonality of the 
explicit Zemike fiinctions which were manually entered. The orthogonality of the explicit Zemike fimctions and the 
Zemike array stmcture was tested next. Finally the orthogonality of the Zemike array stmcture to itself was tested. In 
all three the results were difficult to interpret. The integration is preformed using the qquas routine found in [3], 
p. 148. Because this is a numerical integration we get more and more error as the order increases. Because the first 
part of the test gave similar errors as the second and third, the non-zero terms off of the diagonal may be a result of 
the integration. Better methods of integration will have to be used to validate these conclusions. 

Proving that the transformation algorithms are correct can also be done in an automated manner. To determine 
accuracy of these algorithms we can create simple vectors for each Zemike eigenstate as well as the corresponding 
Taylor amplitude vector. A set of vectors throu^ second order is shown in table 2. We transform the Zemike eigen- 
state into the Taylor basis and compare to the predicted value. To test the reverse transformation, the Taylor ampli- 
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tudes are transformed to the Zemike basis and compared to the expected eigenvector. Since this is all based on integer 
arithmetic, the differences are exactly zero. If there is a non-zero difference then the transformation algorithm has 
been programmed incorrectly. Even more interestingly, by using this method we can determine exactly where the 
problem occurred depending on which transformation jdelded the a non-zero result. 



n m 


Zemike 
Coefficients 


Taylor 
Coefficients 


0 0 


(1) 




1 0 


(0.1,0) 


(0,1,0) 


1 1 


(0,0,1) 


(0,0,1) 


2 0 


(0.0,0,1,0,0) 


(0,0,0.0,2,0) 


2 1 


(0.0,0.0.1.0) 


(-1,0,0,2.0.2) . 


2 2 


(0.0,0.0.0.1) 


(0.0.0.-1.0.1) 



Table 2: Simple Zernike Polynomials and the Equivalent Taylor Monomials 

Testing the precision of these affme transformations can be done by performing two simple tests repeatedly. These 
tests can give a good measure of how accurate these conversions are. To do these tests, create an array of type dou- 
ble to be transformed. The size of the array is (d + l)(c? + 2)/2 where d is the order. Once the array has been created 
with the correct size, fill each element with a random number between -1 and 1 . 

for(int ii = 0; ii < nTerms; ii++) 

initarraytii] = ( (rand( ) / (double)RAND_MAX) - 0.5) * 2;' 

With this array, assume that it contains Zemike amplitudes and perform the transformation to its Taylor equivalent. 
The intermediate array of Taylor amplitudes should then be converted back into its Zemike representation. 

intermediate = fz * initarray; 
finalarray = tz * intermediate; 

In theory, the finalarray should now contain the same values as our initarray. So, by doing this repeatedly 
for each order and with newly created random sequences each time, we can compare the sum of the squares of the dif- 
ferences between initarray and finalarray to determine the accuracy of our transformations. The second test 
is to repeat the above test doing the transformations in reverse. In other words, initarray should be initially 
assumed to be holding the Taylor coefficients. 

intermediate = tz * initarray; 
finalarray = fz * intermediate; 

7. SUMMARY 

The use of circle polynomials of Zemike in optical analysis provides valuable insight. Despite their role in optical 
analysis, computation using the Zemike polynomials is quite inefficient. We outlined how rational polynomials like 
those of Zemike, Legendre, Chebyshev and Laguerre can be handled as afifme combinations of a Taylor monomial 
set. We demonstrated how calculations can be performed much more rapidly in the Taylor basis and how to use inte- 
ger transformations to recover the exact amplitudes in the desired basis. We also discussed some basics of a C++ 
implementation. Finally testing and debugging a large volume of code requires mathematical validation. To insure the 
correctoess of the code, validation was used to ensure the Zemike amplitudes were generated correctly and the trans- 
formations to and from Taylor worked conrectly. Even though Zemike polynomials can tell us a lot about an optical 
system they are not good for doing computation, therefore computation is done in Taylor space. 
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Docket No. WFS.017 
- PATENT - 

[000189] While preferred embodiments are disclosed herein, many variations are possible 
which remain within the concept and scope of the invention. Such variations would become clear 
to one of ordinary skUl in the art after inspection of the specification, drawings and claims herein. 
The invention therefore is not to be restricted except within the spirit and scope of the appended 
claims. 
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Docket No. WFS.017 
-- PATENT - 



CLAIMS 

We claim: 

1. A method of detennining a location of a focal spot on a detector array, comprising: 
defining boundaries of an integration domain for a lenslet; 

performing a center of mass calculation with respect to pixels within the boundaries; and 
correcting for an asymmetry parameter. 
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~ PATENT - 

ABSTRACT OF THE DISCLOSURE 
[000190] A system and method of wavefront sensing with a Shack-Hartmann wavefiront sensor 
precisely locates focal spots o n a detector array, and determines the location of the lenslet array 
with respect to the detector array. 



68 



